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SECTION  I 
INTRODUCTION 

There  is  considerable  current  interest  in  the  effect 
of  nonlinear,  nonconservative  loadings  on  the  stability  of 
elastic  structural  components.  In  particular,  a  new  class 
of  nonlinear,  nonconservative  problems  has  been  introduced 
into  the  field  of  aeroelasticity  with  the  development  of 
efficient  computational  methods  for  unsteady  transonic 
flow.  Nonlinear  transonic  flutter  calculations  pose  a  new 
and  challenging  problem  for  the  aeroelastician  from  both 
an  aeroelastic  analysis  and  an  optimization  standpoint. 

In  this  chapter  some  background  material  is  given  on 
the  nature  of  aeroelastic  instabilities.  This  is  followed 
by  a  statement  of  the  research  objectives.  The  background 
information  points  out  some  of  the  analytical  difficulties 
introduced  into  flutter  calculations  by  nonlinear  transonic 
aerodynamics.  Research  objectives  are  defined  which  are 
designed  to  address  not  only  transonic  nonlinearities,  but 
other  types  of  nonlinearities  as  well.  A  summary  state¬ 
ment  concludes  this  chapter  by  giving  a  concise  statement 
of  the  research  developments  and  applications. 
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1.1  BACKGROUND 


Aeroelasticity  can  be  defined  as  the  study  of  the 
effect  of  aerodynamic  forces  on  a  deformable  body.  In 
general,  these  forces  will  be  nonconservative  and  contain 
both  dissipative  and  circulatory  terms.  Under  certain 
critical  flow  conditions,  the  deformable  body  may  undergo 
large  oscillatory  deformations  which  can  lead  to  structural 
failure  due  to  material  yielding  or  cyclic  fatigue. 

It  is  the  problem  of  the  aeroelastician  to  assure 
that  the  structure  or  structural  component  is  free  from 
failure  when  subjected  to  the  service  flow  conditions.  If 
a  candidate  design  is  not  free  from  failure,  an  appropriate 
design  change  must  be  made  to  provide  for  safety  with 
minimum  loss  in  structural  utility.  The  traditional 
design  changes  are  the  rearrangement  of  structural  stiff¬ 
ness  and  mass  with  total  mass  taken  as  the  measure  of 
structural  utility. 

In  his  classic  text  on  aeroelasticity,  Y.  C.  Fung^ 
distinguishes  between  dynamic  response  problems  and  dynamic 
instability  problems.  The  former  are  characterized  mathe¬ 
matically  as  systems  of  nonhomogeneous  equations  together 
with  prescribed  initial  conditions.  Dynamic  instability 
problems  involving  nonconservative  airloads,  i.e.,  flutter, 
are  characterized  by  systems  of  homogeneous  equations. 
Flutter  instability  occurs  when  an  airflow  parameter  is 


exceeded  as  determined  from  the  solution  of  an  eigenvalue 
problem. 

For  aircraft  flying  in  speed  regimes  other  than 

transonic,  flutter  calculations  are  routinely  formulated, 

The  situation  in  the  transonic  range  is  quite  different; 

however,  since  the  governing  aerodynamic  equations  are 

nonlinear  and  are  quite  difficult  to  solve.  Moving  shock 
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waves  are  very  often  present  which  implies  that  a  non¬ 
linear  mixed  hyperbolic-elliptic  initial  boundary  value 
problem  must  be  solved  to  determine  the  unsteady  airloads. 

One  of  the  most  important  features  of  transonic  flow 
from  the  aeroelastician ' s  viewpoint  is  that  the  airloads 
are  not  related  to  structural  deformation  in  a  linear 
manner.  Dynamic  instability  cannot  be  determined  in 
general  by  solving  an  eigenvalue  problem.  The  problem 
takes  on  some  characteristics  of  a  dynamic  response  calcu¬ 
lation  in  that  initial  conditions  must  now  be  taken  into 
account. 

The  transonic  flow  environment  has  long  presented 
the  designer  with  a  number  of  critical  problems.  Recent 
flutter  model  tests  on  supercritical  wings^  have  shown 
that  catastropic  flutter  occurs  at  the  lowest  flight  speed 
in  the  transonic  Mach  number  range.  It  is  in  this  aero- 
dynamically  nonlinear  range  that  the  structural  designer's 
problem  is  the  most  complex  and  at  the  same  time  the  most 


important.  Lifting  surfaces  and  skin  panels  must  be  made 
free  from  instability  with  minimum  added  weight. 

To  date,  consideration  of  unsteady  aerodynamics  in 
structural  optimization  has  been  limited  to  linear  formula¬ 
tions  of  both  the  structural  and  the  aeroelastic  equations. 

A  number  of  approaches  to  pressure  calculations  based  on 
finite  difference  and  finite  element  solutions  of  the  non¬ 
linear  small-disturbance  velocity-potential  equation  have 
recently  been  developed,  but  flutter  calculations  have  been 
restricted  to  two-dimensional  airfoil  models  having  no  more 
than  three-degrees-of-f reedom.  In  spite  of  the  fact  that 
transonic  flutter  computations  are  still  in  an  early  stage 
of  development,  it  is  appropriate  to  give  consideration  to 
the  establishment  of  optimal  design  methods  for  transonic 
problems . 

All  nonlinear  transonic  flutter  studies  to  date  have 
relied  heavily  on  a  time  history  approach;  i.e.,  numerical 
integration  of  the  equations  of  motion.  This  procedure  has 
been  shown  to  be  quite  time  consuming  and  to  be  subject  to 
numerical  instability.  It  can  be  sometimes  difficult  to 
distinguish  the  numerical  instability  from  a  real  physical 
instability. 

The  method  developed  herein  provides  an  efficient  and 
easy  way  to  implement  an  alternative  to  the  time  history 
procedure.  It  also  has  applicability  to  stability  problems 
other  than  those  induced  by  airflow.  In  addition,  the 


method  is  useful  in  the  optimization  of  system  design 
variables . 

1 . 2  RESEARCH  OBJECTIVES 

The  amplitude  of  disturbance  to  which  a  nonlinear 
aeroelastic  system  is  subjected  can  greatly  affect  its 
stability.  Time  history  solutions  can  be  used  to  assess 
stability,  but  this  approach  can  become  expensive  when 
applied  to  design  optimization  studies.  A  time  history 
solution  would  need  to  be  obtained  at  possibly  a  large 
number  of  steps  in  the  optimization  process.  For  optimi¬ 
zation  purposes,  approximate  methods  from  the  theory  of 
nonlinear  oscillations  are  attractive,  but  their  range  of 
applicability  has  not  been  well  defined  for  use  by  the 
aeroelastician .  Furthermore,  these  methods  have  generally 
been  applied  to  only  single  degree-of-f reedom  systems. 
Clearly,  further  analytical  developments  are  needed  in  this 
area . 

Therefore,  it  is  the  objective  of  this  research  to 
develop  an  efficient  stability  assessment  procedure  for 
systems  subjected  to  nonlinear,  nonconservative  loadings. 
Such  a  stability  assessment  procedure  is  needed  to  solve 
aerodynamically  nonlinear  flutter  and  divergence  problems. 
Although  the  primary  need  for  the  procedure  is  in  tran¬ 
sonic  flutter  studies,  it  is  also  an  objective  of  this 
work  to  apply  the  method  to  other  nonlinear,  nonconsera- 
tive  problems  of  interest  in  analytical  mechanics. 


The  computational  objective  is  to  provide  an  easy 
to  implement  stability  assessment  procedure  which  is  of 
sufficient  efficiency  that  it  can  be  applied  to  nonlinear 
aeroelastic  optimization  problems.  The  end  product  of  this 
research  is  an  approach  to  the  problem  of  determining  the 
effect  of  system  design  variable  changes  on  amplitude  depen¬ 
dent  instability.  To  date,  no  method  has  been  developed 
which  addresses  this  important  design  and  analysis  problem. 

1 . 3  SUMMARY 

The  method  of  imposed  disturbances  is  developed  and 
applied  to  the  problem  of  predicting  amplitude  dependent 
instabilities  in  nonlinear  systems.  For  classical  one- 
degree-of- freedom  systems  such  as  the  van  der  Pol  equation 
or  the  Lewis  servomechanism  equation,  the  load-amplitude 
stability  predictions  agree  very  closely  with  numerical 
integration. 

A  NACA  64A006  airfoil  pitching  in  unsteady  transonic 
flow  is  studied  to  determine  the  system  damping  required  to 
prevent  flutter.  The  trends  predicted  by  numerical  integra¬ 
tion  and  the  imposed  disturbances  method  agree  quite  well 
over  a  wide  range  of  pitch  angle  amplitudes.  The  LTRAN2 
transonic  aerodynamics  code  is  used  to  predict  the  airload 
distributions,  and  it  is  shown  in  the  case  of  increasing 
steady  angles  of  attack,  that  center  of  pressure  shifts  can 
dramatically  affect  divergence.  Airfoil  divergence  occurs 
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when  the  air  speed  dependent  aerodynamic  moment  exceeds  the 
elastic  restoring  moment.  The  critical  air  speed  at  which 
this  occurs  is  the  "divergence  speed,"  and  the  statically 
unstable  airfoil  is  said  to  be  "torsionally  divergent."^ 

Application  of  the  method  of  imposed  disturbances  to 
a  follower  load  problem  is  made  to  demonstrate  the  accuracy 
of  the  method  for  multiple  degree-of-f reedom  systems. 

Again,  the  accuracy  is  quite  good  compared  to  other  methods. 

An  important  static  subcase  of  the  method  is  demon¬ 
strated  in  application  to  a  panel  in  nonlinear  hypersonic 
flow.  A  static  instability  is  predicted  by  the  method  and 
compared  with  Kutta-Merson  numerical  integration.  It  is 
shown  that  the  method  is  over  three  orders  of  magnitude 
faster  than  numerical  integration  and  is  thus  ideally 
suited  for  nonlinear  aeroelastic  optimization  studies  which 
require  a  capability  for  rapid  reanalysis. 
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SECTION  II 


REVIEW  OF  THE  LITERATURE 

Research  in  aeroelastic  optimization  involves 
an  interplay  of  a  number  of  broad  technical  areas.  A 
complete  review  of  the  literature  in  any  single  participat¬ 
ing  discipline  would  be  quite  voluminous.  However,  the 
literature  which  is  relevant  to  the  extension  of  optimiza¬ 
tion  procedures  to  nonlinear  aeroelasticity  is  relatively 
small,  and  it  is  this  subset  which  is  reviewed.  Since  the 
research  needed  for  this  work  grew  out  of  work  in  linear 
aeroelastic  optimization,  the  first  subsection  below  gives 
a  broad  survey  of  work  in  this  area.  This  is  followed  by 
reviews  of  unsteady  transonic  aerodynamics  and  nonlinear 
aeroelastic  stability. 

?.l  AEROELASTIC  OPTIMIZATION 

The  following  survey  is  not  intended  to  provide  an 
exhaustive  review  of  aeroelastic  optimization.  The  excel- 

4 

lent  paper  by  Stroud  gives  a  very  complete  survey  of  work 
done  on  lifting  surfaces.  Some  of  the  more  theoretical 
work  on  panel  flutter  optimization  is  detailed  by  Taylor.^ 
However,  it  is  the  intent  here  to  cite  work  which  has  found 
either  significant  application  or  which  has  provided 


fundamental  theoretical _ insight  into  minimum  weight  design 
with  aeroelastic  constraints.  Reference  is  made  to  the 
paper  by  Lansing,  Lerner  and  Taylor  as  a  detailed  survey 
of  recent  engineering  applications  of  aeroelastic  design 
methods.  The  overall  area  of  structural  optimization  since 
about  1900  is  surveyed  in  References  7  through  10. 

The  term  aeroelastic  optimization  is  used  herein  to 
mean  the  minimum  weight  design  of  a  lifting  surface  or  skin 
panel  subject  to  lower  bound  requirements  on  critical 
flutter  velocity.  Several  other  constraints  are  of  interest 
in  aeroelasticity  such  as  requirements  on  natural  frequen¬ 
cies  and  static  divergence  velocity.  The  review  given 
mainly  deals  with  work  done  on  flutter  constraints,  since 
this  is  the  area  of  prime  concern  in  nonlinear  and  transonic 
aeroelasticity.  Work  relating  to  frequency  constraints  is 
discussed  in  so  far  as  it  can  be  related  to  flutter 
requirements . 

The  early  work  in  aeroelastic  optimization  by 
Turner, MacDonough,  ^  and  Head^^  deals  with  redistribu¬ 
ting  stiffness  to  achieve  a  condition  which  is  indirectly 
related  to  flutter  velocity.  For  example,  the  MacDonough 
condition  is  one  of  uniform  strain  energy  density  under  the 
inertia  loading  associated  with  either  a  wing  bending  or 
torsional  mode.  It  is  necessary,  however,  to  know  before¬ 
hand  the  frequency  requirements  needed  to  prevent  flutter. 
Head'^^  stresses  the  importance  of  the  uniformity 
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of  strain  energy  density  in  the  fundamental  torsional 

mode  as  the  key  to  flutter  prevention  with  minimum  weight. 

14 

The  1968  paper  by  Ashley  and  McIntosh  marks  the 
beginning  of  recent  work  on  aeroelastic  optimization  of 
continuous  systems.  All  known  results  on  optimal  design 
of  continuous  aeroelastic  systems  are  collected  together 
and  it  is  shown  what  functionals  are  uniform  in  minimum 
weight  structures.  Constraints  include  fixed  longitudinal 
and  torsional  frequencies,  torsional  and  chordwise  wing 
divergence  speeds,  torsional  wing  flutter  and  sandwich 
panel  flutter  requirements.  Although  the  sandwich  panel 
flutter  optimization  is  discussed  at  length,  no  numerical 
results  are  presented  in  this  paper. 

The  work  by  Turner^ ^  lays  the  foundation  for  the 
finite  element  approach  to  aeroelastic  optimization.  The 
classical  Lagrange  multiplier  technique  is  used  to  enforce 
the  linear  algebraic  flutter  equation  constraints.  A 
Newton- Raphson  method  is  used  to  satisfy  the  nonlinear 
optimality  equation.  Results  are  given  for  a  three 
element  cantilever  wing  and  for  a  segmented  sandwich 
panel.  The  aerodynamic  theories  used  are  incompressible 
strip  theory  for  the  wing,  and  piston  theory  for  the 
panel  problem.  Numerical  results  for  the  panel  verify 
the  theoretical  prediction  that  the  optimal  panel  thick¬ 
ness  distribution  is  symmetrical  about  its  midpoint. 


The  groundwork  for  continuous  aeroelastic  optimiza-  ' 

tion  laid  by  Ashley  and  McIntosh  is  further  developed  by 
1 6 

Weisshaar.  Solutions  to  one-dimensional  linear  struc¬ 
tural  problems  are  obtained  by  the  transition  matrix 
procedure  adapted  from  the  theory  of  optimal  control. 

For  the  first  time,  quite  detailed  information  is  shown 
on  the  minimum  weight  design  of  a  sandwich  panel  in  super¬ 
sonic  flow.  Weisshaar  also  develops  a  basis  for  discrete 
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element  optimization  via  control  theory  methodology. 

1 8 

McIntosh  gives  quite  a  complete  survey  of  control  theory 

methods  as  applied  to  structural  optimization  problems. 

Optimization  of  discrete  structures  by  numerical 

search  techniques  is  further  extended  by  Rudisill  and 
19 

Bhatia.  The  key  to  their  development  is  the  derivation 
of  expressions  for  partial  derivatives  of  the  flutter 
velocity  with  respect  to  the  member  sizes  in  the  finite 
element  structural  model.  A  search  strategy  is  employed 
based  on  flutter  velocity  gradient  information.  The 
method  is  shown  to  be  quite  computationally  efficient  and 
capable  of  being  used  in  conjunction  with  important 
strength  requirements. 

It  is  interesting  to  note  the  variety  of  design 
sensitivity  methods,  such  as  the  Rudisill-Bhatia  flutter 
derivatives,  that  can  be  used  in  structural  dynamics  and 
aeroelasticity .  In  his  doctoral  dissertation,  van  de 
Vooren^*^  applies  perturbation  theory  to  determine  changes 
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in  flutter  eigenvalues  when  small  changes  in  wing  mass 

and  stiffness  are  made.  Cross  and  Albano'^'*'  further  develop 

the  perturbation  method  and  present  a  working  procedure  for 

rapidly  analyzing  a  large  number  of  wing  store  configura- 

22 

tions  for  flutter  clearance.  Fox  and  Kapoor  compute 
frequency  and  mode  shape  derivatives  for  the  generalized 

algebraic  eigenvalue  problem  with  symmetric  matrices. 

2  2  2  4  2  2 

Rogers,  Rudisill,  and  Rudisill  and  Chu  extend  this 

to  more  general  systems  which  are  not  self-adjoint. 

2  6 

Simpson  studies  the  case  of  derivatives  of  equal  eigen¬ 
values  which  is  applicable  to  the  design  of  systems 
subjected  to  nonconservative  loads  of  the  circulatory 
type.  Both  perturbation  and  derivative  methods  signifi¬ 
cantly  reduce  the  cost  of  reanalysis  when  a  large  number  of 

small  changes  must  be  made  to  a  linear  system. 

27 

Stroud,  Dexter,  and  Stein  develop  a  computer  code 
which  considers  both  flutter  and  strength  in  a  simultane¬ 
ous  manner.  A  plate  model  of  a  delta  wing  is  studied,  and 
design  variables  consist  of  coefficients  in  an  assumed 
polynomial  skin  thickness  distribution.  Move  directions 
for  flutter  resizing  are  generated  by  finite  difference 
approximations  to  flutter  derivatives.  The  overall 
flutter-strength  resizing  is  accomplished  by  an  interior 
penalty  function  mathematical  programming  technique. 

Similar  work  on  multiple  constraints  is  reported  by  Fox, 

2  8 

Muira,  and  Rao  who  treat  flutter,  deflection,  and 
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natural  frequency  constraints  in  a  simultaneous  manner. 

The  work  is  also  notable  since  a  combined  weight-drag 
objective  function  is  minimized. 

29 

In  his  doctoral  dissertation,  Gwin  describes  the 
development  of  a  mathematical  programming  approach  to 
flutter  optimization  based  on  the  method  of  feasible 
directions.  For  the  first  time,  three-dimensional 
aerodynamic  theories  (doublet-lattice  and  Mach-box)  are 
used  to  determine  the  flutter  characteristics.  It  is 
shown  that  flutter  velocity  derivative  information  can  be 
obtained  from  the  rate  of  change  of  complex  eigenvalues 
associated  with  critical  points  on  the  velocity-dam.ping 
diagram..  Using  this  flutter  derivative  information,  Gwin 
efficiently  obtains  minimum  weight  designs  of  several 
complex  finite  element  lifting  surface  models.  Refer¬ 
ence  30  gives  further  details  of  the  results. 

Parallel  to  the  mathematical  programming  methods 

development  are  procedures  based  on  an  energy  related 

optimality  criterion.  One  of  the  most  successful  studies 

31 

is  that  by  Siegel  in  which  the  B-1  wing  structure  is 
resized  to  meet  a  flutter  requirement.  Material  is  added 
to  the  areas  of  maximum  strain  energy  density  in  the  flutter 
mode.  This  is  an  important  advance  over  the  very  early 
flutter  optimization  studies  which  concentrate  on  just  a 
single  free  vibration  mode.  Using  the  Turner  wing  model 
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of  Reference  15,  Bogner  verifies  the  efficiency  and 

utility  of  the  strain  energy  procedure. 

The  flutter  optimality  condition  of  Wilkinson, 

33 

et  al.  requires  uniform  flutter  velocity  derivatives  of 
the  Rudisill-Bhatia  form.  An  extensive  evaluation  of  the 
mathematical  programming  and  energy  based  optimality  pro¬ 
cedures  is  also  presented.  The  result  of  the  effort,  a 
general  purpose  computer  code  FASTOP  (Flutter  And  Strength 
Optimization  Program) ,  capable  of  combining  strength  and 
aeroelastic  optimization,  is  now  available  to  the  aero¬ 
space  community.  A  number  of  significant  detailed  appli¬ 
cations  of  this  code  are  given  in  Reference  34.  It  is 
important  to  note  that  only  linear  aerodynamic  and  struc¬ 
tural  theories  are  included  in  the  code.  The  developments 
reported  herein  are  designed  to  lay  groundwork  for  new 
computer  codes  in  this  area. 

In  order  to  extend  minimum  weight  design  procedures 
to  include  complex  multiple  constraints  and  nonlinear 
effects,  it  is  necessary  to  maintain  and  develop  further  a 
firm  theoretical  base  for  both  the  mathematical  programming 

and  optimality  criteria  approaches.  The  important  work  of 

35 

Prager  and  Taylor  provides  background  for  Plant's  excel- 
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lent  development  on  design  for  specified  critical  loads. 
Reference  36  provides  a  most  complete  mathematical  develop¬ 
ment  of  optimality  for  minimum  weight  design  in  the  presence 
of  circulatory,  gyroscopic,  and  dissipative  forces.  The 


doctoral  dissertations  of  Segenreich,  Johnson,  and  Rizzi^^ 
are  also  prime  examples  of  contributions  which  have  added 
in  a  fundamental  way  to  an  understanding  of  how  mathematical 
programming  and/or  optimality  criteria  can  be  used  in  the 

presence  of  complex  static  and  dynamic  constraints. 

3  8 

In  a  recent  paper  by  McIntosh  and  Ashley,  some 
comparisons  are  made  between  mathematical  programming  and 
optimality  approaches.  Based  on  a  numerical  test  case,  the 
authors  conclude  that  the  optimality  approaches  may  indeed 
be  better  suited  for  incorporation  into  general  design 
procedures  than  certain  mathematical  programming  approaches. 
The  McIntosh  and  Ashley  formulation  allows  for  the  incor¬ 
poration  of  an  "aerodynamic  memory"  effect.  This  formula¬ 
tion  could  be  very  useful  in  the  development  of  a  practical 
aeroelastic  design  capability  which  includes  the  effect  of 
aerodynamic  gust  loadings. 

The  above  survey  shows  that  there  has  been  no  work 
done  to  date  in  transonic  aeroelastic  optimization,  or 
indeed  in  any  area  of  design  optimization  where  the  system 
is  undergoing  nonlinear  motions.  The  next  section  deals 
with  the  status  of  unsteady  transonic  aerodynamics  methods 
which  are  applicable  to  nonlinear  flutter  problems. 

2.2  UNSTEADY  TRANSONIC  AERODYNAMICS 

■j  9 

In  a  recent  survey  of  the  literature,  Borland"^  gives 
an  extensive  review  of  work  relating  to  unsteady  transonic 
flowfield  determination.  A  comprehensive  reference  list 
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shows  that  over  200  papers  and  reports  are  available  in 
the  open  literature  dealing  with  this  subject.  The 
publications  can  be  classified  as  to  solution  method, 
geometry,  time  dependence,  and  type  of  equations  solved. 

Of  the  publications  reviewed  by  Borland,  only  five  are 
listed  in  the  comprehensive  reference  list  as  dealing 
specifically  with  flutter.  All  five  of  the  transonic 
flutter  studies  are  shown  to  be  experimental  in  nature. 

Out  of  a  total  of  171  theoretical  publications,  only  six 
report  data  on  unsteady  pressure  distributions  which  are 
an  essential  ingredient  in  transonic  flutter  calculations. 
Thus,  although  there  is  much  available  literature  on 
"unsteady  transonic  flow,"  only  a  fraction  of  the  litera¬ 
ture  deals  with  methods  directly  applicable  to  aero- 

elastic  calculations. 

40 

McCroskey  makes  the  following  categorization  of 
unsteady  transonic  flow  formulations  in  order  of  increasing 
complexity: 

(1)  linear  potential  theory, 

(2)  local  linearizations, 

(3)  unsteady  linearization  of  small  disturbance 
equations, 

(4)  nonlinear  small  disturbance  equations 

(5)  nonlinear  Euler  equations,  and 

(6)  Navier-Stokes  equations. 


16 


Of  these  formulations,  the  nonlinear  small  disturbance 

methods  have  achieved  a  high  degree  of  development  and 

can  be  used  to  solve  aeroelastic  stability  problems 

while  maintaining  the  nonlinear  nature  of  the  problem. 

It  is  essential  to  use  nonlinear  aerodynamic  theory  to 

look  for  the  onset  of  amplitude  dependent  flutter  and 

2 

divergence.  Tijdeman  shows  essentially  the  same  scheme 
to  categorize  the  analytical  formulations.  It  appears 
that  by  far  the  most  used  computational  method  of  solution 
is  finite  differences.  In  all  six  areas  noted  above, 
much  work  is  reported  in  the  literature  which  makes  use 
of  this  procedure. 

41 

The  monograph  by  Landahl  gives  a  brief  but  straight¬ 
forward  derivation  of  the  basic  transonic  potential  flow 
equation.  The  order  of  magnitude  analysis  needed  to 

arrive  at  the  various  potential  flow  approximations  is 

42  4  3 

presented  by  numerous  investigators  including  Olsen.  ' 

For  an  introduction  to  the  method  of  finite  differences, 

44  4  5 

the  texts  by  Ames  ,  and  Richtmyer  and  Morton  are 

particularly  useful  and  present  much  material  relating  to 

the  solution  of  nonlinear  flow  problems. 

The  alternating-direction  implicit  finite  difference 

46  47 

algorithm  is  used  by  Ballhaus  and  Goorjian  '  to  solve 
the  low-frequency,  small  disturbance  potential  flow 
equations.  A  computer  code,  LTRAN2,  is  now  available 
which  can  predict  the  shock  wave  motions  on  an  oscillating 
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airfoil  or  an  airfoil  with  an  oscillating  flap.  The 

analytical  results  agree  qualitatively  with  other  recent 

48 

experimental  work  of  Tid^eman  and  the  numerical  solutions 

49 

ot  Magnus  and  Yoshihara  who  solved  the  Eulerian  equations 
with  an  explicit  finite  difference  algorithm.  Ballhaus  and 
Goorjian  report  that  a  solution  for  one  cycle  of  unsteady 
motion  can  be  obtained  in  eight  seconds  on  the  CDC  7600 
computer,  compared  to  1500  seconds  for  the  Magnus-Yoshihara 
procedure . 

In  November  of  1978  a  conference  entitled  "Transonic 
Unsteady  Aerodynamics  for  Aeroelastic  Applications"  was 
held  at  Columbus,  Ohio,  under  the  sponsorship  of  the  Air 
Force  Flight  Dynamics  Laboratory.  Papers  presented  at  that 
conference  show  that  much  work  has  been  done  to  predict 
accurately  the  effect  of  a  moving  shock  wave  on  the  pressure 
distribution  at  the  surface  of  an  oscillating  airfoil.  It 
appears  that  three-dimensional  aeroelastic  analyses  will 
soon  be  feasible  on  some  of  the  larger  digital  computers. 

A  full-scale  aeroelastic  optimization  capability  does  not 
appear  possible  in  the  near  future  using  the  analysis  meth¬ 
odologies  that  have  been  proposed  to  date. 

2.3  NONLINEAR  AEROELASTIC  STABILITY 

In  this  subsection  the  area  of  nonlinear  aero- 
elasticity  is  reviewed.  Publications  cited  deal  primarily 
with  either  structural  or  aerodynamic  nonlinearity.  By 
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far  the  most  analytical  work  has  been  done  in  the 
structural  area.  This  type  of  nonlinearity  originates 
usually  in  control  system  linkage  mechanisms  or  from  the 
in-plane  tensile  stresses  in  panels.  The  analytical  work 
that  considers  aerodynamic  nonlinearity  deals  mainly  with 
hypersonic,  transonic,  or  vortex  shedding  phenomena. 

Chapter  10  of  the  text  by  Bisplinghoff  and  Ashley^*^  gives 
an  excellent  introduction  to  some  of  the  more  classical 
methods  used  in  the  analysis  of  nonlinear  systems. 

In  the  paper  by  Woolston,  Runyan,  and  Andrews, 
nonlinear  flutter  calculations  are  performed  on  a  typical 
airfoil  section  which  has  pitch  and  plunge  degrees  of 
freedom.  Structural  nonlinearities  considered  are  due  to 
free-play,  cubic  stiffness,  and  hysteresis  effects 
in  the  restraining  springs.  Numerical  results  from  an 
analog  computer  are  presented  which  indicate  that  stability 
is  highly  dependent  upon  the  magnitude  of  the  initial  dis¬ 
placements.  The  introduction  of  free  play  into  the 
torsional  spring  causes  the  nonlinear  flutter  speed  to  be 
lower  than  the  speed  predicted  without  it.  The  ratio  of 
the  nonlinear  flutter  speed  to  the  linear  speed  decreases 
with  increasing  initial  pitch  displacement.  A  similar 
trend  is  noted  in  the  case  of  a  cubic  variation  in  pitch 
spring  stiffness.  When  hysteresis  is  added  to  the  wing 
control  surface  system,  the  linear  flutter  speed  results 
seem  to  provide  a  good  lower  bound  on  the  nonlinear  case. 
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A  number  of  instances  of  limited  amplitude  flutter  are 

shown.  In  these  situations,  flutter  could  contribute 

significantly  to  fatigue  failure  of  the  system. 

The  "principle  of  harmonic  linearization"  of 

52 

Kryloff  and  Bogoliuboff  is  applied  in  the  paper  by  Shen 

and  Hsu^^  to  the  flutter  problem  of  Reference  51.  By 

balancing  the  fundamental  harmonic,  flutter  velocity 

versus  amplitude  trends  are  shown  which  demonstrate  the 

utility  of  this  procedure.  Woolston  and  Andrews^^ 

point  out  that  direct  integration  of  the  nonlinear  aero- 

elastic  equations  gives  flutter  velocity  as  a  function  of 

the  amplitude  of  the  initial  conditions  of  the  system. 

On  the  other  hand,  the  Kryloff  and  Bogoliuboff  method 

gives  velocity  as  a  function  of  the  limit  cycle  or  long 

55 

term  amplitude  of  the  system.  Shen  elaborates  further 
on  the  harmonic  linearization  method  and  shows  that 
the  long  term  amplitude  results  of  Reference  53  are  not 
directly  comparable  to  those  of  Woolston,  Runyan,  and 
Andrews . 

4 

Transient  and  nonlinear  aeroelastic  effects  are 
considered  by  Brunelle^^,  and  an  approximate  method  is 
developed  to  handle  systems  with  one  or  two  degrees  of 
freedom.  Brunelle  points  out  that  the  Kryloff- 
Bogoliuboff  solution  assumes  that  only  weak  nonlinearity 
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and  small  damping  are  present.  He  then  sh.ows  how  to 
solve  the  van  der  Pol  equation  starting  with  an 
assumed  solution  which  contains  an  exponential  time 
function  to  account  for  strong  nonlinearities.  In 
the  two-degree-of-f reedom  case,  the  Brunelle  method 
becomes  quite  lengthy.  Results  for  the  van  der  Pol 
equation  agree  quite  well  with  numerical  integration 
and  indicate  that  this  approximate  method  can  predict 
accurately  the  limit  cycle  behavior  of  the  can  der  Pol 
equation  with  strong  nonlinearity.  No  numerical  results 
were  given  for  the  two-degree-of-freedom  system. 

Nonlinear  aeroelastic  response  of  panels  is  discussed 
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at  length  by  Dowell  in  Chapter  3  of  his  monograph.  The 
basic  technique  of  modal  expansion  is  presented  and 
applied  to  skin  panels  where  in-plane  stretching  is  the 
source  of  the  nonlinearity.  The  effects  of  static 
pressure  differential,  in-plane  loading,  and  curvature 
on  flutter  dynamic  pressure  are  presented.  In  References 
58  and  59,  Dowell  also  discusses  direct  numerical  inte¬ 
gration  of  the  nonlinear  panel  flutter  equations  from  a 
given  set  of  initial  conditions. 

Perturbation  and  harmonic  balance  methods  are  applied 

6  0 

by  Kuo,  Morino,  and  Dugundji  to  large  deflection  panel 
flutter  problems.  Results  agree  well  with  those  obtained 


by  the  direct  integration  method.  It  is  shown  that  struc¬ 
tural  damping  can  have  a  destabilizing  effect  on  flutter 
at  high  mass-air  ratios.  Both  the  perturbation  method 
and  the  harmonic  balance  method  are  shown  to  require  only 
about  four  percent  of  the  computing  time  needed  for 
numerical  integration. 

6 1 

The  direct  method  of  Liapunov  has  been  studied  by 

62  63 

a  number  of  investigators  including  Parks,  George,  and 
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Dimantha  and  Roorda.  As  applied  to  panel  flutter.  Parks 
shows  that  the  Liapunov  method  can  lead  to  highly  conserva¬ 
tive  bounds  on  critical  dynamic  pressure.  Such  results 
can  be  highly  dependent  on  the  form  of  the  Liapunov 
functional  that  is  constructed.  The  work  of  Diamantha  and 
Roorda  deals  with  this  problem,  and  they  show  how  to 
construct  functionals  for  some  simple  nonconservative 
systems.  The  power  of  the  Liapunov  method  lies  in  the 
fact  that  stability  "in  the  large"  can  be  determined. 

The  problem  of  constructing  appropriate  functionals  for 
complex  nonlinear  systems  such  as  those  encountered  in 
aeroelasticity  remains  an  area  of  current  research. 

On  a  practical  engineering  level,  the  effect  of 
structural  nonlinearities  on  flutter  behavior  is  of  con¬ 
tinued  interest  to  the  aircraft  designer.  The  recent  papers 
65  66 

by  Breitbach,  Laurenson  and  j.rn,  and  Sensburg  and 
Schoen^^  show  that  there  is  continuing  concern  within  the 


aerospace  industry  over  effects  of  free  play  and  nonlinear 
stiffness  on  missile  and  aircraft  control  surface 
flutter  boundaries.  Approximate  analytical  tools  are 
being  built  around  the  method  of  harmonic  balance,  and 
experimental  models  are  being  used  to  check  the  accura  . 
of  the  analyses. 

Work  in  the  area  of  aerodynamically  nonlinear  flutter 

appears  to  be  quite  limited.  Eastep  and  Mclntosh^^ ' 

investigate  the  stability  of  a  flat  panel  in  nonlinear 

hypersonic  flow.  Both  numerical  integration  and  the 
64 

Dimantha-Roorda  energy  considerations  lead  to  the  con¬ 
clusion  that  a  panel  can  become  unstable  in  hypersonic 
flow  if  the  energy  input  of  the  initial  disturbance  is 
sufficiently  large. 

Instabilities  of  a  nonlinear  aerodynamic  nature  can 
also  involve  separated  flow  around  an  oscillating  airfoil. 
At  relatively  low  subsonic  speeds,  a  complex  interaction 
involving  the  fluctuating  lift  due  to  vortex  shedding  and 
an  elastic  structure  can  occur.  This  "vortex-induced 
vibration"  is  self -excited  and  self-limited.  A  precise 
modeling  of  the  flow  environment  would  require  solving 
the  full  Navier-Stokes  equations,  but  some  investigators 
such  as  Hartlen  and  Currie^^  have  proposed  an  attractive 
alternative.  In  their  study,  use  is  made  of  experimental 
lift  data  for  a  long  circular  cylinder  in  low-speed  flow 
to  derive  a  nonlinear  lift-oscillator  model  of  the 
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aerodynamic  environment.  The  governing  equations  for 

the  system  can  be  shown  to  be  equivalent  to  the  van  der 

Pol  oscillator.  Approximate  solutions  are  obtained  based 

on  the  approximate  method  given  by  Minorsky.”^^  Agreement 

with  experiment  is  quite  good;  however,  other  investi- 

40 

gators  such  as  McCroskey  have  criticized  the  lift- 

oscillator  model  on  the  grounds  that  it  is  not  derived 

from  the  fundamental  principles  of  fluid  mechanics. 

A  related  aerodynamically  nonlinear  problem,  stall 

7  2 

flutter,  is  investigated  by  Dugundji  and  Chopra.  In  this 

primarily  experimental  study,  the  stability  of  two- 

dimensional  wings  is  investigated  at  high  angles  of  attack 

in  subsonic  flow.  Nonlinear  divergence  is  shown  to  occur 

at  angles  of  attack  above  about  10  degrees.  In  the  case 

of  bending-torsion  flutter,  limit  cycle  behavior  is  shown 

at  torsional  amplitudes  as  great  as  40  degrees.  The  onset 

of  this  coupled  mode  flutter  can  be  predicted  with 

reasonable  accuracy  using  linear  aerodynamic  theory. 

In  the  transonic  speed  range,  a  single  degree-of- 

freedom  instability  known  as  "aileron  buzz"  has  been  known 

to  occur.  This  phenomenon  results  from  a  separated  flow 

condition  and  is  very  difficult  to  predict  analytically. 

As  early  as  1946,  these  unstable  aileron  oscillations  were 

7  3 

of  concern  to  the  aeroelastician .  Smilg  presents  a 
practical  engineering  approach  to  the  problem  and  reports 
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in  detail  on  the  effect  that  various  airfoil  and  aileron 

parameters  have  on  the  oscillations.  As  a  result  of  this 

74 

investigation  and  also  a  later  paper  by  Smilg,  the 
importance  of  single  degree-of-freedom  aeroelastic 
stability  considerations  is  well  established. 

Within  the  last  two  years,  several  important  tran¬ 
sonic  flutter  studies  have  been  completed  which  incorporate 

state-of-the-art  unsteady  transonic  aerodynamic  methods. 

47 

Ballhaus  and  Goor^ian  present  flutter  results  of  a 

NACA  64A006  airfoil  pitching  in  transonic  flow.  The 
46 

LTRAN2  code  is  coupled  with  a  simple  numerical  integra¬ 
tion  procedure  to  solve  the  structural  and  aerodynamic 

75-77 

equations.  Rizzetta  also  investigates  the 

use  of  numerical  integration.  The  equations  of  motion  for 
a  three -degree-of-freedom  airfoil-aileron  system  are 
coupled  with  LTRAN2  aerodynamic  lift  and  moment  calcula¬ 
tions  and  integrated  simultaneously.  Rizzetta  shows  that 

7  8 

the  Adams-Moulton  implicit  predictor-corrector  method 
has  superior  numerical  stability  compared  to  other  pro¬ 
cedures  tried,  such  as  the  Euler -Cauchy  and  the  Milne 
three-point  implicit  schemes.  A  number  of  amplitude 
dependent  flutter  cases  are  shown  which  indicate  the 
essential  nonlinear  character  of  the  transonic  flutter 
problem . 
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7  9  8  0 

The  work  by  Yang,  Striz,  and  Guruswamy  '  repre¬ 
sents  a  comprehensive  study  of  the  effect  of  various 
structural  and  aerodynamic  parameters  on  the  flutter 

behavior  of  two-dimensional  airfoils  in  transonic  flow. 

8 1 

The  results  of  Reference  79  indicate  that  relaxation  and 
47 

indicial  aerodynamic  methods  give  comparable  flutter 

solutions  in  the  low  transonic  Mach  number  range.  In 

Reference  80,  a  direct  integration  method  is  coupled  with 

the  LTRAN2  code,  to  compute  time-history  flatter  solutions. 

Results  for  a  NACA  64A006  airfoil  agree  quite  well  with 

47 

those  of  Ballhaus  and  Goorjian.  Time-history  solutions 
and  the  indicial  and  relaxation  methods  give  significantly 
different  estimates  of  aerodynamic  coefficients.  Refer¬ 
ence  80  states  that  the  differences  are  due  to  the  time 
linearization  assumption  of  the  indicial  and  relaxation 
methods . 

Qualitative  results  on  the  influence  of  shock  wave 

motion  on  typical-section  wing  flutter  are  presented  by 
8  2 

Ashley.  Shock  parameters  considered  are  displacement 
amplitude,  phase  lag  behind  quasi-steady  behavior,  steady 
strength,  and  chordwise  position.  The  effect  of  the  shock 
is  shown  to  be  generally  destabilizing  for  single-degree-of- 
freedom  pitching  motion.  For  bending-torsion  cases  the 
shock  can  be  either  stabilizing  or  destabilizing,  depending 
on  the  system  parameters.  Qualitative  results  such  as 
Ashley's  should  provide  insight  into  the  basic  transonic 


flutter  mechanism  which  will  be  of  great  use  in 
further  analytical  work. 
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SECTION  III 


THEORETICAL  DEVELOPMENT 

The  "method  of  imposed  disturbances"  is  derived  based 

p  O 

on  the  concept  of  conservation  of  energy  in  a  limit  cycle. 

For  certain  problems,  nonlinear  instability  occurs  at  zero 
frequency  and  the  method  takes  the  form  of  a  set  of  gener¬ 
alized  force  balance  equations.  The  limitations  of  the 
Krylof f-Bogoliubof f  method  are  shown  in  the  first  subsection 
below.  This  is  followed  by  a  general  formulation  of  the 
method  of  imposed  disturbances.  It  is  shown  how  the 
method  is  easily  used  in  conjunction  with  sophisticated 
aerodynamic  calculations.  As  an  alternative  to  the  cumber¬ 
some  time  history  approach  to  nonlinear  flutter  calculations, 
the  method  is  highly  efficient  and  gives  clear  physical 
insight  into  the  nature  of  the  instability. 

3.1  LIMITATIONS  OF  THE  KRYLOFF-BOGOLIUBOFF  METHOD 

The  literature  survey  shows  that  the  harmonic  balance 
method  due  to  Kryloff  and  Bogoliuboff  (K-B  method)  is 
currently  being  applied  to  obtain  approximate  solutions  to 
nonlinear  flutter  problems.  In  what  follows,  two  simple 
examples  will  be  studied  to  point  out  the  utility  and  also 
the  limitations  of  this  approach.  The  discussion  of  the 
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method  follows  that  of  Minorsky®^  and  the  original  Kryloff- 

52 

Bogoliuboff  reference. 

3.1.1  Basic  Considerations 

Consider  a  differential  equation  for  a  single 
degree-of-f reedom  dynamic  system  which  has  the  form 


(3.1) 


X  +  2^(jj^x  +  ujj^x  =  ijf(x,x) 


where  f  is  in  general  a  nonlinear  function  of  the  displace¬ 
ment,  X,  and  the  velocity,  x.  The  quantities  and  c, 
are  the  undamped  natural  frequency  and  damping  factor, 
respectively.  The  nonlinear  terms  are  scaled  by  the  small 
parameter  y.  Equation  (3.1)  can  be  written  as  a  pair  of 
first  order  equations 


(3.2a) 


^  =  -F (x,x) 


(3.2b) 


where 


(3.3) 


F(x,x)  =  +  (iJ^X  -  uf(x,x) 


Division  of  (3.2a)  by  (3.2b)  gives 


(3.4) 


Q(x,y) 

P(x,y) 


where  y  = 


(3.5) 


Q(x,y)  =  -F(x,y) 


and , 


(3.6) 


P(x,y)  =  y 


If  X  and  y  are  considered  to  be  points  in  a  plane,  i.e., 
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the  phase  plane,  then  it  can  be  shown  that  Equation  (3.4) 

has  only  one  integral  curve  passing  through  each  regular 

point  of  the  plane.  A  point  in  the  phase  plane  is  regular  if 

P  and  Q  do  not  both  vanish;  otherwise,  the  point  is  said  to 
34 

be  singular.  A  limit  cycle,  C,  is  a  special  solution  of 
Equation  (3.2)  with  the  following  properties;  (1)  it  is  a 
closed  trajectory,  and  (2)  all  other  solutions,  C',  either 
approach  C  as  t-^<»  or  as  t-»-  -<»,  if  c'  approaches  C  as  t-*-®, 
then  C  is  a  stable  limit  cycle,  and  in  the  case  that  C 
approaches  C  as  t-^  the  limit  cycle  is  said  to  be  unstable. 
Not  all  nonlinear  systems  exhibit  limit  cycle  behavior; 
however,  the  phenomenon  does  occur  in  many  instances  which 
are  of  practical  concern  in  analytical  mechanics.  Typical 
limit  cycle  behavior  in  the  phase  plane  is  shown  in  Figure  1. 
Stable  limit  cycle  behavior  is  exhibited  in  Figure  la,  since 
both  representative  trajectories  converge  toward  the  same 
closed  curve.  Just  the  opposite  is  true  in  Figure  lb  as  both 
curves  diverge  away  from  the  closed  path,  C.  Trajectories 
which  start  exactly  on  the  unstable  limit  cycle  will  continue 
indefinitely  on  this  path.  The  unstable  limit  cycles  are 
difficult  to  calculate  numerically  since  a  small  accumulation 
of  error  will  deflect  the  trajectory  either  to  the  origin  or 
to  some  unbounded  path. 


(x'x)jrl  =  x'rn  +  x'^n!?^ 


3.1.2  Expressions  for  Amplitude  Derivatives 

The  first  approximation  method  of  Kryloff- 
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Bogoliuboff  gives  information  concerning  the  ^  =  0  solu¬ 
tion  of  Equation  (3.1).  A  solution  is  assumed  to  be  of  the 
form 

(3.7)  X  =  asin(cot  +  <})) 

where  a  and  4)  are  slowly  varying  functions  of  time  and  w 
is  constant  and  equal  to 

The  stability  of  motion  is  closely  related  to  the 
value  of  the  time  derivative  of  the  amplitude.  If  a  is 
a  monotonically  decreasing  function  of  time,  then  the 
motion  is  stable.  Instability  occurs  when  a  is  always 
positive.  Explicit  expressions  for  a  can  give  much  insight 
into  the  stability  of  the  system. 

A  requirement  that  the  velocity  have  the  form 

(3.8)  X  =  awcos(cot  +  <p} 
gives  the  condition  from  Equation  (3.7)  that 

(3.9)  asin((jjt  +  41)+  a^cos  (wt  +  4)  =  0 

When  Equation  (3.8)  is  differentiated  and  substituted  into 

Equation  (3.1),  a  second  condition  on  a  and  is  obtained. 

•  *  • 

(3.10)  acocos((ot  +  4)  “  aai4sin(cjt  +  4)  =  ijf(x,x) 

where  the  arguments  of  f  are  given  by  Equations  (3.7)  and 
(3.8)  . 

Equations  (3.9)  and  (3.10)  are  easily  solved 

1 


to  get 


(3.11a) 


a  =  uf(x,x)  cos(cjt  +  (J)/a) 


(3.11b)  $  =  -Mf(x,x)  sin(a)t  +  !p)/aw 

Since  the  right  hand  sides  of  Equation  (3.11)  are  almost 
periodic  functions,  they  can  be  represented  in  the  first 

approximation  as  a  Fourier  series  of  the  form 

00 

(3.12a)  Uf(x,x)  COSY  =  I 

^  n=l 

+  P '  (a)  sin  ny] 
n 

00 

(3.12b)  uf(x,x)  siny  =  ^ 

+  Q^(a)  sin  ny] 

where  y  =  a)t  +  (}> .  Substitution  of  Equations  (3.12)  into 
(3.11)  gives  expressions  which  can  be  integrated  over  the 
time  interval  t  to  t  +  T  where  T  is  the  period.  Considering 
a  and  $  to  be  constant  whtin  this  interval  leads  to 
(3.13a)  [a(t+T)  -  a(t)]/T  =  uPQ(a)/aj 

(3.13b)  [(i>(t+T)  -  (t)(t)]/T  =  uQQ(a)/aw 

If  the  period  is  great  enough,  the  difference  quotients  on 
the  left  sides  of  Equation  (3.13)  can  be  replaced  by  time 
derivatives.  Thus  using  the  usual  Fourier  expressions  for 
Pq  and  Qq ,  the  first  approximation  of  the  K-B  theory  is 
obtained  as  follows 
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f2lT 

>  (3.14a)  a  =  (u/2-n<jj)  |  f(a  sinY,aa)  cosy)  cosydy 

^  o 

.  (-2^ 

4i  =  -y/27Tcoa)  f(a  sinY,aa)  cosy)  sinydy 

o 

where  a  and  d  are  treated  as  constants  in  the  integration. 

^  The  power  of  this  method  can  be  seen  in  an 

I  application  to  the  van  der  Pol  equation 

► 

[  (3.15)  X  +  X  =  y (1-x^) X 

F*  In  this  case  the  time  rate  of  change  of  amplitude  as  compu¬ 

ted  by  Equation  (3.14a)  is 

(3.16)  a  =  (ua/2)  (l-a^/4) 

It  is  suggested  by  this  equation  that  for  positive  y,  the 
motion  tends  toward  a  stable  limit  cycle  of  amplitude  2.0 

(see  Figure  la).  For  negative  y,  a  is  negative  at  ampli- 

♦ 

tudes  less  than  2.0  and  positive  for  those  greater  than  2.0. 
Hence  in  this  case  an  unstable  limit  cycle  is  indicated  (see 
Figure  lb) .  Such  conclusions  for  the  van  der  Pol  equation 
are  verified  by  numerical  integration  in  Chapter  4.  The  K-B 
theory  is  shown  to  give  accurate  trend  information  concern¬ 
ing  limit  cycle  behavior  for  even  moderate  values  of  y  (about 
1.0)  . 

The  method  just  outlined  gives  equivalent 
results  to  the  variation  of  parameter  method  described  by 
Cunningham®^  and  the  van  der  Pol  method  (see  Stoker,®® 
pp.  149-153) .  In  the  latter  method,  a  solution  of  the  form 
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(3.17) 


X  =  a^(t)  sincat  +  a2(t)  cosojt 

is  assumed.  Such  a  form  can  be  obtained  from  Equation 

(3.7)  by  elementary  trigonometry. 

3.1.3  The  Harmonic  Balance  Form  of  the  K-B  Method 

The  harmonic  balance  form  of  the  K-B  method 
54 

used  by  Shen  in  his  nonlinear  flutter  work,  takes  on  a 
slightly  different  but  equivalent  form.  Again  using  the 
^  =  0  form  of  Equation  (3.1),  the  first  step  is  to  assume 
a  series  solution 

00 

(3.18)  X  =  ^  a  sin  nwt 

n=l 

and  substitute  it  into  the  given  equation.  This  leads  to 
00 

(3.19)  y  (-w^n^  +  b  )  sin  nwt  =  0 

n=l 

where  the  b  constants  are  the  Fourier  sine  coefficients 
n 

of  (jo^x  =  yf(x,x).  In  explicit  form,  the  first  coefficient 
is 

rlTr/ui  2 

(3.20)  ^1  ~  ^  -  uf(x,x)]  sinwt  dt 

*  o 

where  the  first  term  of  the  series  Equation  (3.18)  is  used 
to  compute  the  integral.  Balancing  the  first  harmonic  in 
Equation  (3.19)  gives 

(3.21)  u)^  =  b^/a^ 

where  b^  =  b^(a^)  is  computed  from  Equation  (3.20).  The 
first  order  solution  is  thus 
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(3.22) 


mm 


(3.22)  X  =  sin [b^ (a^) t/a^] 

3.1.4  The  K-B  Estimate  of  Stability  for  the  Reduced 
Pendulum  Equation 

As  an  example  of  harmonic  balance  procedure, 
consider  the  reduced  nonlinear  pendulum  equation  (see 
Reference  52,  p.  14) 


(3.23) 


X  +  (g/l)  (x-x'^/6)  =  0 


where  x  is  the  angular  displacement  of  the  pendulum  and  g 
and  i  are  the  acceleration  due  to  gravity  and  the  pendulum 
length,  respectively.  Equation  (3.23)  gives  accurate 
results  up  to  displacement  angles  of  about  30  degrees.  It 
will  be  shown  below  that  larger  angles  can  result  in  an 
instability.  Although  unstable  motion  is  not  observed  in 
the  physical  system,  it  serves  to  point  out  a  difficulty 
with  the  unjudicious  application  of  the  harmonic  balance 


method . 


A  straightforward  application  of  Equations 


(3.18)  through  (3.20)  yields  the  following  amplitude  depen¬ 
dent  frequency  equation 


(3.24) 


=  (g/£)  (l-a^/8) 


From  this  equation,  it  is  concluded  that  real  frequencies 
can  be  computed  for  all  amplitudes  of  motion  up  to  a^  =  /8. 
Beyond  that  point.  Equation  (3.24)  has  complex  roots  and  the 
displacement  solution  of  Equation  (3.22)  has  a  term  which 
grows  exponentially  with  time. 


Figure  2  shows  the  solution  of  Equation  (3.23) 

8  7 

that  is  obtained  by  the  Kutta-Merson  integration  technique . 

A  slightly  more  general  equation  has  been  included  in  the 
computer  code,  but  the  case  shown  in  Figure  2  is  for 
?  =  =  0 ,  p  =  g/Si  =  1.0  with  an  initial  displacement  of 

X  =  2.448.  The  motion  is  stable  and  has  aperiodic  frequency 
of  0.2527  radians  per  second  which  compares  very  well  with 
the  value  of  0.2509  radians  per  second  obtained  from  Equa¬ 
tion  (3.24)  with  a^^  =  2.448.  At  a  slightly  higher  initial 
Hisplacement  of  2.449,  the  motion  is  divergent  as  shown  in 
Figure  3.  Thus  Equation  (3.24)  overestimates  the  critical 
initial  displacement  which  causes  instability  by  about 
15  percent.  It  is  also  important  to  note  that  the  K-B 
estimate  of  rate  of  change  of  amplitude  using  Equation 

(3.14a)  gives  a  E  0  for  all  amplitude  levels. 

The  K-B  method  can  be  applied  to  an  n-degree- 

of-freedom  system;  however,  the  expressions  for  the  rate  of 

change  of  amplitude  require  n  multiple  integrations  of  the 

88 

nonlinear  terms.  Such  a  calculation  for  prediction  of  aero- 
dynamically  nonlinear  flutter  could  be  quite  formidable  when 
explicit  expressions  for  the  aerodynamics  are  not  available. 

An  additional  difficulty  with  the  K-B  method  is 
that  strong  linear  damping  terms  cannot  be  included  on  the 
left  hand  side  of  the  equations.  As  stated  in  the  litera¬ 
ture  survey,  Brunelle^^  modified  the  K-B  method  to  include 
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Figure  2.  Stable  Nonlinear  Motion  of  Reduced  Nonlinear 
Pendulum  Equation  (Kutta-Merson  Integration) 


Figure  3.  Unstable  Nonlinear  Motion  of  Reduced  Nonlinear 


Pendulum  Equation  (Kutta-Merson  Integration) . 


strong  damping  terms  and  nonlinearities;  however,  the  exten¬ 
sion  of  the  method  to  multiple  degrees-of-f reedom  systems 
appears  difficult. 

3.2  THE  METHOD  OF  IMPOSED  DISTURBANCES 

In  the  last  section  it  was  shown  how  an  approximate 
method  due  to  Krylof f-Bogoliubof f  could  be  used  to  gain 
insight  into  the  limit  cycle  behavior  of  systems  such  as 
those  governed  by  the  van  der  Pol  equation.  For  other 
problems  such  as  the  reduced  nonlinear  pendulum,  the 
method  can  lead  to  quite  unconservative  results.  This 
section  describes  an  approximate  method  which  can  deal  with 
systems  in  which  the  strict  application  of  the  K-B  method 
can  lead  to  inaccuracies.  The  method  proposed  herein  is  an 
improvement  over  the  usual  K-B  method  in  that  it  is  now 
possible  to  account  for  strong  linear  damping  and  nonlinear 
static  instabilities. 

Before  starting  the  development,  several  properties 
of  nonlinear,  nonconservative  systems  are  reviewed.  This 
is  followed  by  the  statement  of  the  design  optimization 
problem  for  nonlinear  systems  and  the  approach  to  its 
solution . 

3.2.1  The  Load-Amplitude  Diagram 

The  primary  concept  to  be  reviewed  is  that  of 
the  load-amplitude  diagram.  A  mechanical  system  governed 
by  the  following  set  of  matrix  equations  is  to  be  studied: 
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(3.25)  Mx+Cx-pGx+Kx=  /f(x,x) 

where 

M  =  generalized  mass  matrix 
C  =  generalized  damping  matrix 
G  =  generalized  geometric  stiffness  matrix 
K  =  generalized  stiffness  matrix 
X  =  generalized  coordinates 
f  =  load  vector  which  is  nonlinear  in  the 
generalized  coordinates  and  velocities 
p  =  static  load  parameter 
A  =  dynamic  load  parameter 

There  may  often  exist  a  critical  value  of  the  dynamic  load 
parameter,  A  =  A,  beyond  which  the  motion  of  the  system 
increases  without  bound-  A  further  important  feature  of  the 
system  is  that  A  can  depend  upon  the  magnitude  of  the 
initial  disturbance.  It  is  sometimes  possible  to  define  a 
single  disturbance  parameter.  A,  and  to  numerically  solve 
Equation  (3.25)  at  a  given  A.  For  example,  in  nonlinear 
panel  flutter  problems  the  quantity  A  can  be  taken  to  be 
the  maximum  amplitude  at  the  three-quarters  chord  point. 

Thus  (A, A (A))  plots  can  be  made  as  shown  in  Figure  4  by 
varying  A  until  self-sustained  motion  is  achieved. 

Figure  4a  depicts  what  can  be  termed  a  "hard 
instability"  since  increasing  the  amplitude  of  the  distur¬ 
bance  also  increases  the  critical  value  of  the  load.  Thus 
there  is  an  apparent  increase  in  stiffness  in  the  system. 
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Figure  4.  Possible  Load-Amplitude  Relationships  for  a  Nonlinear 
Nonconservative  System. 


system  of  Equations  (3.25),  deleting  the  nonlinear  terms 
contained  in  the  load  vector.  At  some  value  of  greater 
than  Xq,  i.e.,  as  shown  in  Figure  4a,  the  system  is 
unstable  at  low  amplitudes  of  disturbance.  At  higher  dis¬ 
turbances  the  system  is  again  stable.  In  the  case  of  a 
single  degree-of-f reedom  system,  the  phase  plane  representa¬ 
tion  would  show  the  existence  of  stable  limit  cycle  motion. 
From  an  energy  standpoint,  at  low  amplitude  the  system  is 
gaining  energy  from  the  dynamic  loads,  but  at  higher  ampli¬ 
tudes  the  system  is  losing  energy,  and  the  motion  cannot  be 
maintained . 

A  soft  instability  is  shown  in  Figure  4b.  In 
this  case,  a  linearized  analysis  gives  a  critical  value  of  X 
which  is  an  upper  bound  on  the  stable  values  of  the  loading 
parameter.  At  some  value  of  X  less  than  Xq,  the  system  is 
unstable  until  a  critical  amplitude  is  reached  and  then 
unstable  motion  occurs.  In  the  phase  plane,  unstable  limit 
cycle  motion  would  be  observed.  It  is  this  case  of  soft 
instability  which  is  of  great  technological  importance.  A 
linearized  analysis  proves  to  be  highly  misleading  if  large 
amplitude  disturbances  are  expected  to  be  applied  to  the 
system.  Figure  4c  shows  a  possible  behavior  which  is  a 
combination  of  both  situations  just  described. 

3.2.2  Optimization  of  Load-Amplitude  Curves 

At  this  point  the  (A, A)  design  optimization 
problem  is  discussed.  As  seen  from  Figure  4b,  designing 


the  system  so  that  is  greater  than  some  fixed  value  is 
inadequate.  It  is  apparent  that  the  stability  of  the  system 
must  be  assured  within  some  rectangular  domain  bounded  by 
both  loading  and  amplitude  requirements.  Figure  5  shows 
possible  "optimized"  load-amplitude  plots  for  these  cases. 

By  changing  system  design  variables  in  the 

generalized  matrices,  it  is  possible  to  reconfigure  the 

load-amplitude  curves.  In  Figure  5,  the  required  load,  a„, 

R 

and  the  required  amplitude,  A  ,  block  off  a  region  of  the 
diagram  in  which  it  is  desired  to  have  only  stable  combina¬ 
tions  of  A  and  A.  The  situation  shown  in  Figure  5a 
dictates  that  system  changes  be  made  which  primarily  affect 
the  A  requirement.  In  practice,  it  would  not  be 
immediately  obvious  what  changes  to  make;  however,  the 
following  comments  concerning  what  is  called  here  "design 
in  the  A-mode"  are  in  order. 

A  first  order  approximation  to  the  solution  of 

the  design  problem  of  Figure  5a  is  entirely  covered  within 

the  framework  of  methodology  already  developed  for  the 

optimization  of  linear  systems.  For  example,  if  Figure  5a 

is  the  load-amplitude  plot  for  a  panel  in  two-dimensional 

nonlinear  flow,  then  A  is  the  dynamic  pressure  parameter, 

and  A  can  be  taken  as  the  panel  deflection  amplitude  at  a 

selected  point.  In  nonlinear  panel  flutter  work,  plots 

5  8 

like  this  are  frequently  developed.  Thus,  in  order 
to  achieve  the  X_  requirement,  the  equations  for  the 

R 
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Figure  5.  Optimized  Load-Amplitude  Relationships 


panel  are  first  linearized,  and  then  structural  resizing 
performed  based  on  a  standard  mathematical  programming 
method,  for  example.  Reference  5  contains  a  survey  of  the 
methods  that  have  been  used  to  optimize  linear  panels  to 
achieve  a  flutter  requirement. 

Such  an  approach  is  not  without  pitfalls.  It 

is  very  possible  that  the  optimized  panel  would  then  have 

the  amplitude  problem  shown  in  Figure  5b.  This  situation 

must  be  treated  as  a  separate  design  problem.  The  potential 

difficulty  with  a  A  to  A  mode  switch  is  similar  to  the 

mode  switching  problem  encountered  in  linear  aeroelastic 

optimization.  The  work  by  O'Connell,  Hassig  and 
8  9 

Radovcich  discusses  switching  in  flutter  modes  and  the 
occurrence  of  "hump  modes." 

What  appears  to  be  lacking  completely  is  an 
approach  to  reconfigure  the  A-mode  instability  shown  in 
Figure  5b.  The  more  complex  nonlinear  "hump  modes"  of 
Figure  5c  can  be  treated  as  a  combination  of  A-mode  and 
A-mode  optimizations.  It  is  therefore  the  specific  objec¬ 
tive  of  this  research  to  develop  an  approximate  method  to 
determine  sensitivity  of  points  on  the  A-mode  instability 
plot  to  changes  in  the  design  variables.  When  this  informa¬ 
tion  is  available,  it  is  then  entirely  within  the  state  of 
the  art  to  develop,  for  example,  a  gradient  based  optimiza¬ 
tion  procedure.  The  objective  is  to  make  the  smallest 


possible  changes  in  the  system  design  variables  which 

remove  the  A(A)  curve  from  the  required  instability  region, 

R  (shaded  area  in  Figure  5)  . 

The  development  of  a  method  to  compute  A-mode 

sensitivities  is  the  nonlinear  analog  to  the  key  development 
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in  the  dissertation  by  Gwin.  In  the  case  of  linear 
flutter,  he  derives  sensitivity  information  for  selected 
points  on  a  V-g  (velocity-damping)  diagram.  Moves  are  made 
in  design  space,  based  on  these  sensitivities,  which  recon¬ 
figure  a  specific  flutter  mode  so  that  the  system  is  stable 
below  a  required  velocity.  With  the  sensitivity  calculation 
in  hand,  the  most  difficult  part  of  the  optimization  problem 
is  solved. 

3.2.3  Mathematical  Development  of  Sensitivities 

Sensitivities  for  the  A-mode  optimization  are 
now  defined  in  more  mathematical  terms.  It  is  assumed  that 
there  exists  a  set  of  independent  design  variables,  h, 
associated  with  the  system  described  by  Equations  (3.25). 
When  a  specific  set  of  these  variables  is  known,  the 
generalized  mass,  damping,  and  stiffness  matrices  can  be 
determined:  i.e., 

M  =  M(h)  ,  C  =  C  (h)  ,  K  =  K(h) 

The  response  of  the  system  is  a  function  of  h,X,  the  initial 
conditions,  and  time. 


u^,  I  .  .Ill, 


X  =  xCh,A,XQ,XQ;t) 

Let  (Xq,Xq)  denote  a  specific  set  of  initial  conditions  which 
results  in  a  self-sustained,  periodic  motion  of  the  system. 

So  at  steady -motion, 

(3.26)  X  =  x(h,XQ,XQ;t) 

(3.27)  A  =  A(S,Xq,Xq) 

(3.23)  oj  =  n(h,XQ,XQ) 

(3.29)  T  =  T(h,Xg,Xg)  =  2tt/Q 

where  w  and  T  are  frequency  and  period,  respectively. 

The  characteristic  amplitude.  A,  is  defined  as 

(3.30)  A  =  max  [X^ (h,XQ ,X  ; t) ] 

i ,  t£l 

where  X^  are  the  components  of  X  and  I  is  the  time  interval 
for  one  period  of  motion.  Thus  X  has  the  form 

(3.31)  X  =  A^ 

where  1'  is  the  "normalized"  coordinate  vector  at  the 
steady-state  condition.  In  general, 

(3.32)  ?  =  (h,XQ,XQ;t) 

and 

(3.33)  A  =  A(h,XQ,XQ) 

The  required  sensitivities  of  points  on  the  A- 
mode  instability  boundary  can  be  represented  symbolically  as 
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(3.34) 


A,^  =  G (h,XQ,XQ)  i  =  1, . . . ,Np 


where  (  ),^  denotes  partial  differentiation  with  respect  to 
t  h 

the  i^  component  of  h,  and  N  is  the  total  number  of  design 

P 

variables.  Equation  (3.34)  can  be  obtained  by  differentia¬ 
tion  of  Equation  (3.27);  however,  that  functional  relation¬ 
ship  is  unknown. 

It  will  now  be  shown  how  an  exact  expression 

for  the  right  hand  side  of  Equation  (3.27)  can  be  derived. 

i-T 

Premultiply  the  governing  Equations  (3.25)  by  X  and  inte¬ 
grate  over  one  cycle  of  motion  to  get  at  the  critical 
condition 


(3.35) 


Pi  i-T  1  -*-T  -*•  1  -►T  Pi 

X^M  X  -  I  pX^  G  X  +  i  X"  K  Xj  = 

I  •  •  •  * 

f  (A  X  r {X,X)  -  X  C  X)  dt 


The  fact  that  M,  G,  and  K  are  symmetric  has  been  used  to 
compute  the  left  side  of  Equation  (3.35).  Define  the  non¬ 
conservative  work  done  during  the  steady-state  motion  as  the 
right  side  of  Equation ( 3 . 35 ) , 


(3.36; 


NC  r  Pt  -►I 

r  =  1__AX  f(X,X)  -  X  C  xJ  dt 

J  n 


Note  that  the  conservative  terms  on  the  left  side  of 
Equation  (3.35)  are  identically  equal  to  zero.  This  is  a 
result  of  the  fact  that  conservative  terms  evaluated  at  the 
upper  time  limit  have  the  same  value  as  at  the  lower  limit. 
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tude  of  steady-state  oscillation.  This  assumption  does  not 
limit  the  applicability  of  the  method.  It  just  provides  an 
unambiguous  way  to  define  amplitude  (see  Equation  (3.30)). 
From  the  study  of  limit  cycle  behavior,  it  is  felt  that  it 
is  reasonable  to  effectively  remove  the  initial  conditions 
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from  the  problem,  since  the  long  term  amplitude  (bounded  or 
unbounded)  can  be  achieved  by  any  numiser  of  initial  con¬ 
ditions.  Hence  Xq  and  are  not  single  valued  functions 
of  A.  It  is  the  character  of  the  normalized  response,  1', 
that  is  important.  By  treating  A  as  a  scale  factor  on  , 

the  relation  A  =  A(A,?,S)  can  be  obtained  from  the  condition 
NC  ... 

W  =0.  This  is  valid  if  Assumption  2  holds.  Again,  by 
studying  typical  nonlinear  systems,  it  can  be  found  that 
the  character  of  the  normalized  response  and  period  do  not 
greatly  change  over  a  wide  range  of  system  parameters. 

Hence,  Assumption  3  is  justified.  Assumption  4  is  valid  as 
long  as  the  motion  is  steady-state.  If  the  period  becomes 
infinite,  a  static  subcase  must  be  considered.  This  is 
presented  in  Section  3.2.6. 

It  is  very  important  to  keep  in  mind  that  a 
rapid  analysis  tool  is  being  developed  here  for  application 
to  optimization  studies.  A  completely  rigorous  analysis 
step  is  not  intended  through  the  entire  resizing  process. 

It  is  felt  that  the  assumptions  used  herein  are  reasonable 
considering  the  purpose  for  which  the  method  is  intended. 

3.2.5  Statement  of  the  Method 

The  method  of  imposed  disturbances  can  now  be 
stated.  For  a  nonlinear,  nonconservative  system  under¬ 
going  steady-state  oscillations,  the  work  done  on  the 
system  by  the  nonlinear  external  forces  must  be  balanced 
by  the  work  done  by  the  internal  damping  of  the  system. 
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This  is  represented  by  the  condition  =  0.  Given  a 

normalized  response  mode,  'J',  which  is  a  kinematically 
admissible  motion,  the  amplitude  parameters.  A,  can  be 
varied  independently  to  achieve  =  0.  This  expression 

can  then  be  solved  for  the  critical  load,  A.  The  first 


result  for  steady-state  motion  is 


(T  •  '  (T 

(3.38)  A(A,S)  =  A  l'^C(h)f  dt/  1''^  l(A,?,?)  dt 

0  0 


and  the  sensitivities  are 


(3.39)  A,^(A,h)  =  A 


dt/  7'^  f(A,?,?)  dt 


The  conclusion  is  that  damping  changes  are  the  main  quan¬ 
tities  to  be  varied  in  order  to  remove  steady-state  A-mode 
instabilities . 

No  application  of  a  procedure  equivalent  to 
the  method  of  imposed  disturbances  has  been  found  in  the 
literature.  Numerical  studies  are  shown  in  Chapter  4  which 
indicate  that  the  assumptions  of  Section  3.2.4  are  not  too 


severe.  The  method  fills  the  important  technological  gap 
that  exists  between  methods  that  are  used  for  nonlinear 
analysis  and  those  that  are  used  for  optimization. 

3.2.6  The  Static  Solution  Subcase 

An  important  subcase  can  be  found  when  the 
nonlinear  external  forces  have  the  form 


(3.40)  f(X,X)  =  ?^(X)  +  f2(X,X) 
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Let  ^  =  '{'^  be  an  unknown  constant  vector.  When  this  is 
substituted  into  Equation  (3.25)  together  with  Equa¬ 
tion  (3.40),  the  following  static  force  balance  equations 
are  obtained: 

(3.41)  F  H  -  [K(h)  -  pG(h)]  =  0 

For  some  given  A  =  A^  and  A  =  7^  0,  a  solution  of  these 

nonlinear  algebraic  equations  gives  the  static  equilibrium 
solution  't'^.  Thus  the  entire  A(A,S)  can  be  generated  by 
solving  Equation  (3.41)  in  a  parametric  manner.  Since 
stiffness  terms  appear  in  Equation  (3.41),  it  can  be  con¬ 
cluded  that  changing  these  terms  is  important  in  removing 
static  A-mode  instabilities. 

3.2.7  Use  of  the  Method 

In  the  next  chapter  a  series  of  numerical 
examples  are  presented  which  show  that  the  steady-state  and 
static  A  (A)  curves  computed  by  the  method  of  imposed  dis¬ 
turbances  agree  quite  well  with  numerical  integration.  The 
method  of  imposed  disturbances  is  particularly  useful  for 
aerodynamically  nonlinear  flutter  and  divergence  calcula¬ 
tions.  In  these  cases,  explicit  forms  for  ^  ana  are  not 
always  known,  but  may  be  evaluated  instead  by  numerically 
solving  some  complex  flow  equation.  By  imposing  a  steady- 
state  disturbance  on  the  wing  panel,  the  calculation  of  the 
W*  involves  only  integrations  and  matrix  multiplications. 
For  flutter  calculations,  no  interaction  with  the  structural 
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part  of  the  equations  is  necessary  if  'i'  is  properly 
selected.  The  aerodynamic  calculations  are  performed 
using  ?  to  determine  the  imposed  time  dependent  flow 
boundary  conditions.  The  next  chapter  shows  that  good 
results  can  be  obtained  when  is  computed  by  just  solving 
the  linearized  problem. 

The  procedure  for  flutter  calculation  is  thus 
to  solve  either  a  low  amplitude  nonlinear  problem,  or,  if 
possible,  to  solve  a  fully  linearized  problem.  Then  this 
solution  is  used  as  an  imposed  disturbance  as  input  to  a 
series  of  flow  calculations.  Each  flow  calculation  involves 
a  different  scale  factor.  A,  on  the  mode  4^.  The  approximate 
flutter  solution  is  at  that  amplitude  which  causes  the  net 
nonconservative  work  to  vanish.  The  method  of  imposed 
disturbances  is  proposed  as  an  analysis  tool  which  is 
complementary  to  time  history  flutter  solutions.  It  can 
also  give  a  first  order  estimate  of  the  sensitivity  of  the 
system  to  various  design  changes. 
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SECTION  IV 
NUMERICAL  RESULTS 

In  this  chapter  the  method  of  imposed  disturbances 

is  applied  to  a  series  of  nonlinear  problems.  The  first 

section  below  describes  the  application  of  the  method  to 

several  widely  studied  nonlinear  equations.  Comparisons 

are  made  with  various  numerical  integration  methods. 

Nonlinear  flutter  and  divergence  calculations  are  then 

studied  which  use  the  transonic  aerodynamic  code,  LTRAN2 , 

to  compute  the  air  loadings.  The  chapter  is  concluded  with 

a  follower  load  problem  and  a  hypersonic  panel  flutter 

problem  to  demonstrate  the  applicability  of  the  method  to 

multiple  degree-of-f reedom  systems. 

4.1  NONLINEAR  OSCILLATOR  STUDIES 

Two  widely  studied  nonlinear  differential  equations 

are  the  van  der  Pol  equation  and  the  Lewis  servomechanism 

equation.  Accurate  numerical  solutions  are  given,  for 

90 

example,  in  the  dissertation  by  Knese.  In  what  follows, 
the  method  of  imposed  disturbances  is  shown  to  give  trend 
results  which  agree  well  with  parametric  studies  performed 
using  numerical  integration. 
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4.1.1  The  van  der  Pol  Equation 

The  usual  form  of  the  van  der  Pol  equation  does 
not  contain  a  strong  linear  damping  term.  It  will,  however, 
be  included  here  to  demonstrate  the  applicability  of  the 
method  of  imposed  disturbances  to  such  problems.  Thus 
including  the  damping  term  the  equation  is 

(4.1)  X  +  2(;a)^x  +  =  y(l-x^)x 

The  notation  is  the  same  as  that  used  in  Section  3.1. 

An  imposed  disturbance  is  given  to  the  system 

of  the  form 

(4.2)  X  =  A4' 

where 

(4.3)  4'  =  sinf^t 


This  is  used  in  the  total  nonconservative  work  expression, 


(4.4) 


|^M(1-X^)X  -  2ccoj^x]x  dt 


which  can  easily  be  computed  as 

NC  2  2 

(4.5)  =  (TTufiA^/4)  [(8w^/u)  (u/2oj^-C)-A^] 

It  is  observed  that  if  u>0  and  c<u/2a}^,  then 
NC  2 

W  <0  whenever  A  >  (Sajj^/y)  (u/2co^-C)  ;  i.e.,  the  motion  is 

2 

stable.  If  the  same  condition  holds  on  c,  but  A  is  less 
than  ( 8a)j^/y )  (y/2wj^-? )  then  the  motion  is  unstable.  The 
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optimum  value  of  ^  which  will  preclude  instability  at  any 
amplitude  is 

(4.6)  i;  =  u/2co^  {y>0) 

and  in  which  case  the  motion  is  always  of  decreasing 
amplitude . 

Figure  6  shows  phase  plane  results  for  the  damped 
van  der  Pol  equation  (NDET=1)  using  the  Kutta-Merson  numeri¬ 
cal  integration  method  (NAM=1) .  The  step  size  (H)  is  0.05 
with  500  steps  (NSTEP)  taken.  The  initial  displacement 
(Xl(0))  is  0.01  and  the  initial  velocity  (X2(0))  is  zero. 

Value  of  u(XMU),  C(ZETA),  and  (OMEGA)  are  1.0,  0.0,  and 
1.0,  respectively.  In  this  case.  Equation  (4.5)  predicts 
W  >0  until  an  amplitude  of  2.0  is  reached.  As  shown  in 
Figure  6,  numerical  integration  indicates  that  the  displace¬ 
ment  increases  initially  and  is  limited: to  about  2.0  in 
absolute  value  for  large  time.  , 

In  Figure  7,  an  initial  displacement  amplitude 
of  3.3  is  given  and  the  motion  decreases  and  again  reaches 

a  limiting  amplitude  of  2.0.  Such  behavior  is  consistent 

NC 

with  the  observation  that  W  <0  for  this  condition. 

To  further  check  the  validity  of  Equation  (4.5), 
assume  that  it  is  required  that  sufficient  damping  be  added 
to  the  system  which  will  make  it  stable  for  all  (u,A)  in 
the  range  -0.2<u<0  and  0<A</^.  Solving  Equation  (4.5)  for 
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NDET  -  1  NAM 

-  1  H  -  .  05 

ic's  XI C05.  X2ca5 

«  .  01  0.  00 

XMU.  ZETA.  OMEGA  - 

I.  00  0.  0000  1.  00 

CPTIME  -  .46 

NS7EP  -  500 

Figure  6.  Stable  Limit  Cycle  Behavior  of  van  der  Pol 
Equation  (Small  Disturbance) . 
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XDOT 


NDET  -  1  NAM  -  i  H  -  .  05 
IC'S  XI  C05.  X2<05  -  3.30  0.00 

XMU,  ZETA.  QMEGA  -  1.  00  0.  0000  1.  00 

CPTIME  -  .  47  NSTEP  -  500 


Figure  7.  Stable  Limit  Cycle  Behavior  of  van  der  Pol 
Equation  (Large  Disturbance)  . 
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NC 

Z  at  the  W  =0  condition  gives 

(4.6)  ^  ~  y  (4-A^) /8(jo^  (y<0) 

so  for  the  conditions  indicated,  c  must  be  no  smaller  than 
0.05.  Figure  8  shows  the  phase  plot  for  an  initial  dis¬ 
turbance  just  slightly  less  than  The  motion  is  indeed 

stable  for  =  0.05.  However,  in  Figure  9  the  initial 
disturbance  is  just  slightly  more  than  and  the  motion  is 
unstable  at  the  same  Thus  the  accuracy  of  Equation  (4.6) 

for  this  case  is  verified. 

4.1.2  A  Modified  van  der  Pol  Equation  for  Load- 
Amplitude  Solution 

Consider  a  modified  form  of  the  van  der  Pol 

equation 

(4.7)  X  +  +  oj^x  =  A(l+x^)x 

where  the  scale  factor,  X,  on  the  nonlinear  terms  is 
considered  to  be  a  dynamic  loading  parameter. 

For  vanishing  nonlinearity,  the  critical  load  is 

(4.8)  ^ 

and  a  solution  is 


(4.9) 


=  A  t 
o  o 
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Figure  8.  Unstable  Limit  Cycle  Behavior  of  Damped 

van  der  Pol  Equation  (Small  Disturbance) . 


NOET  -  1  NAM  -  1  H  -  .  05 
IC'  S  XI  C0>  ,  X2  <05  -  2.  50  0.  00 

XMU.  ZETA.  OMEGA  -  20  .  0500  1.  00 

CPTIME  -  .17  NSTE?  -  500 


Figure  9.  Unstable  Limit  Cycle  Behavior  of  Damped 

van  der  Pol  Equation  (Large  Disturbance) . 
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where 


(4.10)  4*^  =  sintD^t 

Equation  (4.9)  is  the  zero  amplitude  solution  of  Equation 
(4.7),  and  hence  an  imposed  disturbance  of  the  form  of 
Equation  (4.2)  can  be  used  to  compute  the  total  nonconserva 
tive  work.  This  leads  to 

(4.11)  =  (A^tt/J^)  [A(1+A^/4)  - 

NC 

At  the  condition  of  neutral  stability,  W  =  0,  an 
appropriate  relation  between  critical  load  and  amplitude  is 

(4.12)  A  =  2coj^/(1+A^/4) 

Alternatively,  A  in  terms  of  A  is 

(4.13)  A  =  2  /(2Cajj^-A)/A 

NC 

Curves  of  W  =0  are  shown  on  Figure  10  for 
various  values  of  Time  history  integration  results  are 

shown  to  verify  the  accuracy  of  Equation  (4.12) . ,  All 
"K-M"  points  are  the  result  of  Kutta-Merson  numerical  inte¬ 
gration  using  x(0)  =  A  and  x(0)  =  0  as  the  initial  condition 
Equation  (4.13)  predicts  the  shape  of  the  load-amplitude 
curve  quite  accurately  and  gives  the  proper  trend  as  ?  is 


increased. 


4.1.3  The  Lewis  Servomechanism  Equation 

In  control  theory,  a  nonlinear  equation 

referred  to  as  the  Lewis  Servomechanism  equation  is 
90 

encountered.  For  the  purposes  herein,  it  is  useful  to 
write  it  in  the  form 

•  2  >1 

(4.14)  X  +  2<;(jj^x  +  ca^x  =  X{]xl  -  l)x 

Again  using  Equation  (4.2)  as  the  initial  disturbance,  the 
total  nonconservative  work  is 
Nr  2 

(4.15)  W  (A  J^A/2)  [4A/3-tt(1+2cw^/A] 

NC 

For  small  A,  W  is  negative  and  for  large  A 
it  is  positive.  Hence  the  existence  of  an  unstable  limit 
cycle  is  indicated.  This  behavior  is  verified  by  Kutta- 
Merson  numerical  integration  of  the  equation  as  shown  in 
Figures  11  and  12.  In  Figure  11,  the  phase  plane  plot  for 
the  case  A  =  1.0,  C  =  0,  =  1.41  is  shown  where  a  pure 

initial  velocity  condition  of  3.30  is  given  the  system. 
These  conditions  produce  stable  motion  of  decreasing  arrpli- 
tude;  at  a  velocity  condition  of  3.40  the  system  is  unstable 
as  shown  in  Figure  12.  Only  the  initial  velocity  condition 
is  different  between  these  two  plots. 

At  the  neutrally  stable  condition  Equation 
(4.14)  can  be  solved  for  A  in  terms  of  A,!ii^,  and  ;  to  get 

(4.16)  A  =  (3TT/4)  (l  +  2^Wj^/A) 
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NDET  -  3  NAM 

-  1  H  -  .  05 

IC'S  XI  C0?.  X2C05 

K  0.  00  3.  30 

XMU.  ZETA.  CMEGA  - 

1.00  0.0000  1.41 

CPTIME^-  .46 

NSTE.5^  -  500 

11.  Unstable  Limit  Cycle  Behavior  of  Lewis 

Servomechanism  Equation  (Small  Disturbance) • 


NDET  -  3  NAM 

-  1  H  " 

.  Z5 

IC'S  XI  CZ5,  X2C0? 

V  z.  zz 

3.  40 

XMU.  ZETA.  OMEGA  - 

1.  zz  z. 

ZZZ0  1.  41 

CPTIME  -  .  IS 

NSTEP  - 

5ZZ 

Figure  12.  Unstable  Limit  Cycle  Behavior  of  Lewis 

Servomechanism  Equation  (Large  Disturbance) 
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In  a  manner  identical  with  the  modified  van  der  Pol 


equation  study,  load-amplitude  plots  for  various  values 
of  ;;  can  be  computed  by  Equation  (4.16)  and  compared  to 
numerical  integration.  These  comparisons  are  shown  in 
Figure  13.  As  with  the  previous  study,  the  accuracy  of 
the  method  of  imposed  disturbances  is  quite  good  for  a 
wide  range  of  parameters. 

4.2  TRANSONIC  DIVERGENCE 

Torsional  divergence  of  a  wing  is  a  static 
instability  that  occurs  when  the  destabilizing  aerodynamic 
moment  exceeds  the  stabilizing  elastic  restoring  moment. 
Figure  14  shows  a  typical  section  airfoil  with  an  aero¬ 
dynamic  moment  acting  on  it.  The  aerodynamic  moment 
can  be  expressed  as 


(4.17) 


1  2  2 
M,,  =  V  c  C 
0  2  oo  <»  m 


where  P  is  the  free-stream  air  density,  V  is  the  free- 

stream  velocity,  c  is  the  airfoil  chord,  and  C  is  the 

m 

nondimensional  aerodynamic  moment  coefficient.  In  low- 
speed  flow  the  aerodynamic  moment  coefficient  can  be 
accurately  represented  as  a  linear  function  in  angle  of 
attack ,  Uq ;  i .e . , 


(4.18) 


m 


C  a 


ma 


0 
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4  6  8  10  12  14  K 

Load-Amplitude  Relationship  for  Lewis 


Servomechanism  Eauation. 


Figure  14.  A  Two-Dimensional  Airfoil  in  Steady  Flow 


If  Oq  ia  too  large  or  if  the  flow  is  transonic,  this 
linear  relation  may  no  longer  hold  over  the  range  of 
of  interest. 

The  onset  of  divergence  occurs  when 

(4.19)  M-  =  M 

Q  e 

where  is  the  elastic  restoring  moment  which  can  be 
expressed  as 

(4.20)  M  =  K  ct,, 

e  a  0 

It  is  assumed  that  the  torsional  spring  is  linear;  hence 

the  spring  constant,  K  ,  is  not  a  function  of  a,. 

'a  0 

Equation  (4.19)  can  b-?  rewritten  as 

(4.21)  =  C 

K  m 

where  is  the  nondimensional  spring  coefficient 

(4.22) 

The  LTRAN2  transonic  aerodynamic  code  was  used  to 
investigate  the  case  when  C  is  a  function  of  ctn  •  Steady 
pressure  coefficient  distributions  about  a  NACA  64A006 
airfoil  were  calculated  at  a  free  stream  Mach  number  of 
0.88.  A  finite  difference  grid  containing  113  chordwise 
points  and  97  vertical  points  was  used  to  solve  the  flow 
problem. 
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Table  1  shows  the  LTRAN2  lift  coefficients  and 
centers  of  pressure  that  were  computed  for  various  angles 
of  attack.  The  lift  coefficient  is  defined  as 

(4.23)  -  P^) 

where  P  =  P_  (x)  and  P  =  P  (x)  are  the  pressures  on  the 

Xi  L>  U  U 

lower  and  upper  surface,  respectively,  of  the  airfoil. 

The  moment  coefficient  can  be  computed  from 

(4.24) 

where  E,.  and  C  are  the  locations  of  the  spring  pitch  axis 

J  o 

and  center  of  pressure,  respectively,  expressed  as  a 
fraction  of  the  chord. 

Figures  15  through  21  show  the  steady  pressure 
coefficient  distribution,  Cp ,  on  the  upper  and  lower 
surfaces  for  angles  of  attack  between  zero  and  one-half 
degrees.  The  critical  pressure  coefficient  is  also 
shown . 

If  the  spring  pitch  axis  is  assumed  to  be  located  at 
.586,  it  is  observed  that  at  about  =  .25°  the  center  of 
pressure  is  located  at  the  pitch  axis.  For  smaller  angles 

it  is  forward  of  this  axis  and  for  larger  angles  it  is  to 


TABLE  1 


LIFT  COEFFICIENTS  AND  CENTERS  OF  PRESSURE 
FOR  NACA  64006  AIRFOIL  IN  STEADY  FLOW 

MACH  NUMBER  =  .88 


a 

^L 

.05 

.080743 

.567475 

.10 

.160636 

.567611 

.15 

.232955 

.576185 

.20 

.286957 

.582854 

.25 

.325688 

.586350 

.  30 

.358254 

.589188 

.35 

.386346 

.591300 

.40 

.408527 

.592814 

.45 

.430277 

.593120 

.50 

.447601 

.593464 

.75 

.501767 

.586323 

1.00 

.533986 

.567874 

1.50 

,578334 

.  540518 

Pitch  axis,  X/C  =  .586 


Figure  16.  Steady  Pressure  Coefficient  Distribution  64A006  Airfoil 


Figure  18.  Steady  Pressure  Coefficient  Distribution,  64A006  Airfoil 


Figure  19,  Steady  Pressure  Coefficient  Distribution,  64A006  Airfoil, 


steady  Pressure  Coefficient  Distribution  64A006  Airfoil, 
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Figure  21.  Steady  Pressure  Coefficient  Distribution,  64A006  Airfoil 


the  aft.  Although  this  pitch  axis  location  may  not  be  very  j 

I 

realistic,  this  example  serves  to  highlight  a  key 
mechanism  of  nonlinear  divergence;  namely,  that  center  of 
pressure  location  is  a  function  of  angle  of  attack.  This 
phenomenon  is  not  observed  in  linear  flow  problems. 

The  shift  in  center  of  pressure  is  directly  due  to 
the  pressure  of  shock  waves  on  the  upper  and  lower  surfaces. 

These  pressure  discontinuities  move  toward  the  trailing  edge 
on  the  upper  surface  as  angle  of  attack  increases,  and 
toward  the  leading  edge  on  the  lower  surface.  This  produces 
a  lateral  change  in  the  pressure  distribution  which  grows 
faster  than  the  transverse  changes  which  are  occurring  toward 
the  leading  edge.  It  is  this  feature  of  pressure  distribu¬ 
tion  which  causes  the  moment  coefficient  to  be  nonlinear. 

The  fact  that  C  is  nonlinear  can  be  seen  from  Equation 

(4.24)  since  both  center  of  pressure  and  lift  are  functions  ^ 

of  angle  of  attack  in  this  case. 

Figure  22  shows  a  graphical  solution  to  Equation 
(4.21).  A  total  of  three  quilibrium  positions  are  observed: 

(1)  the  origin,  (2)  =  +.20°,  and  (3)  =  -.20°.  The 

origin  is  an  unstable  equilibrium  position,  since  for  very  j 

small  Uq  the  destabilizing  aerodynamic  moment  is  larger  in 
absolute  value  than  the  restoring  moment  due  to  the  spring. 

For  positive  such  as  A^  on  Figure  22,  the  motion  tends  to 

I 

( 


I 

t 


increase  toward  the  stable  equilibrium  position  A.  At 
positions  larger  than  about  .20®,  the  motion  decreases 
toward  A.  At  negative  a^,  poiint  B  is  a  stable  position  and 
motions  started  at  points  and  B2  will  quickly  stabilize 
toward  B . 

In  this  divergence  study  the  imposed  disturbance  is  a 
constant  and  stability  is  determined  from  a  simple  moment 
balance  relation.  This  is  a  single  degree-of-f reedom 
example  of  the  static  solution  subcase  discussed  in 
Section  3.2.6.  However,  moments  instead  of  forces  are 
balanced  for  the  transonic  divergence  problem. 

4.3  SINGLE  DEGREE-OF-FREEDOM  TRANSONIC  FLUTTER 

Time  history  flutter  analysis  results  are  obtained 
for  an  airfoil  pitching  in  transonic  flow.  Results  obtained 
by  numerical  integration  are  compared  to  those  of  the 
method  of  imposed  disturbances.  The  optimal  amount  of 
system  damping  needed  to  prevent  flutter  is  computed  by 
both  methods  as  a  function  of  amplitude.  The  numerical 
integration  and  imposed  disturbance  methods  both  predict 


similar  damping-amplitude  trends. 

4.3.1  System  Geometry  and  Equation  of  Motion 


Figure  23  shows  the  geometry  of  a  NACA  64A006 


pitching  in  transonic  flow.  The  streamwise  coordinate, 
is  nondimensionalized  by  the  chord  c.  A  pitch  axis  loca¬ 
tion  of  C  =  0.5  has  been  selected  for  this  study. 


The  equation  of  motion  is 
(4.25)  la  +  ga  +  K^a  =  M 

where  I,  g,  and  are  the  system  pitch  inertia,  damping, 

and  spring  stiffness,  respectively,  and  M  is  the  unsteady 

aerodynamic  moment.  Following  the  nondimens ionalization  of 

4  7 

Ballhaus  and  Goorjian,  the  Equation  (4.25)  can  be  rewritten 
as 


(4.26) 

where 


a"  +  A,a'  +  A_a  =  A.C 
1  2  3  m 


(4.27a)  A^  =  g/cjl 

2 

(4.27b)  A-  =  K  /(D  I 

2  a 

(4.27c)  A3  = 

and  (  )  '  denotes  differentiation  with  respect  to  cot  where  (o 

is  an  arbitrary  frequency  that  has  been  introduced  to  make 

the  time  scale  nondimens ional . 

The  transonic  moment  coefficient,  C  ,  is  a 

m 

complicated  function  of  Mach  number,  cot,  a,  a',  and  reduced 
frequency,  k^.  The  reduced  frequency  definition  is  based 
on  chord;  i.e.. 
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set  up  to  compute  unsteady  lift  and  moment  coefficients 

for  an  imposed  harmonic  airfoil  motion.  To  compute  a  time 

history  solution  in  which  the  airfoil  is  free  to  interact 

with  the  structural  term  in  the  equation,  it  is  necessary 

to  add  a  numerical  integration  algorithm  to  the  code. 

91 

Goorjian  uses  the  following  explicit  finite 

difference  scheme  to  solve  the  same  pitch  flutter  prob- 

lem*^^  that  is  presented  here.  Assume  that  at  the  n^^ 

time  step,  a,  a',  and  C  are  known  and  have  values  denoted 

m 

by  cx  ,  a',  and  C  .  The  expression  for  the  second 
n  n  mn 

derivative  of  a  at  this  time  step  is  computed  by  solving 
Equation  (4.26)  to  get 


{4.29a) 


a"  =  -A, a'  -  A_a  +  A-C 
n  In  2  n  3  mn 


The  following  two-term  Taylor  series  expansion  for  a  gives 


in  terms  of  known  quantities 


(4.29b) 


a  ,  ,  =  a  +  ha'  +  -^h^a" 
n+1  n  n  2  n 


where  h  is  the  time  step.  The  one-term  Taylor  series 
for  a ' 


(4.29c) 


a '  +  ha " 

n  I 


is  the  final  difference  equation  that  is  needed.  The 
LTRAN2  code  computes  after  and  have 

been  determined  from  Equations  (4.29). 


4.3.3  Transonic  Pitch  Flutter  Analysis  by  the  Method 
of  Imposed  Disturbances 
Let 

(4.30)  a  =  a^siny 


be  an  imposed  airfoil  motion  where  y  =  wt  and  is  the 
amplitude  of  pitching  oscillation.  The  total  nonconservative 
work  done  by  the  aerodynamic  and  dissipative  terms  during 
one  period  of  motion  is 


(4.31) 


W 


,NC 


C2tt 

[A-C  (a,a' ,y) -A, a ']  a'  dy 
j  Q  j  m  j. 


where  the  imposed  disturbance  of  Equation  (4.30)  is  used  to 
compute  and  the  other  a  related  terms  in  the  integral . 

Using  Equation  (4.30)  in  Equation  (4.31)  and 
rearranging  terms  gives 


(4.32) 


ttNC, 

W  /oL^-rr  =  a^A^  - 


where  a^^  is  the  usual  Fourier  cosine  coefficient 


(4.33) 


^1  TT 


Yj_+2Tr 


C  (a ,a '  ,y)  cosydy 
m 


The  time  value,  y  is  taken  to  be  6 it  in  the 
numerical  calculations.  This  is  required  in  order  to 
establish  the  periodicity  of  the  airfoil  pressure  distribu¬ 
tion.  The  stability  criterion  is  thus 

>0  unstable 

(4.34)  ^  =0  neutrally  stable 

<0  stable 
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NC 

Solving  W  =0  for  gives  the  required  damping  needed 
to  prevent  flutter  at  an  imposed  disturbance  amplitude  of 
Thus  the  required  damping  is 

(4.35)  Aj^  =  a^^  A^/a^ 

In  the  next  subsection  specific  data  will  be  used  to  compare 
the  results  of  this  equation  with  numerical  integration. 

4.3.4  Numerical  Results  for  Pitching  Airfoil 

All  the  following  calculations  have  been  per¬ 
formed  at  a  Mach  number  of  0.88  and  a  reduced  frequency  of 
0.1.  These  appear  to  be  typical  values^ which  can  lead 
to  nonlinear  flutter  behavior  at  moderate  angles  of  attack. 

Figure  24  shows  the  time  history  results  for  C 

m 

and  a  using  the  LTRAN2  aerodynamic  code  for  the  case  of 

=  .25°,  A^  =  .99,  A2  =  1.2278,  A^  =  0.9596.  Only  and 
A^  are  variables  in  this  study.  All  cases  have  been  run  for 
the  values  of  A2  and  A^  shown  above.  Four  complete  cycles 
of  imposed  motion  are  given  the  airfoil  in  order  to  establish 
a  periodic  airflow  condition.  This  is  followed  by  approxi¬ 
mately  two  cycles'of  free  motion  as  computed  by  the  finite 
difference  method  of  Section  4.3.2.  The  initial  conditions 


for  the  free  motion  are  the  a,  a',  and  at  the  end  of  the 
fourth  cycle  of  imposed  motion.  Figure  24a  shows  the  complete 
time  history  and  Figure  24b  shows  a  full  scale  plot  of  the 
free  motion.  Figure  25  shows  the  time  history  of  the  lift 


Imposed 


Figure  24a.  Moment  and  Angle  of  Attack  Time  History  for  NACA  64A006  Airfoil 
M  =  .88,  a,  =  .25°,  A  =1.05,  A„  =  1.2275,  A^  =  0.9596. 


£5000  *  00E 


Figure  24b.  Moment  and  Angle  of  Attack  Time  History  for  NACA  64A006  Airfoil 

=  .88,  =  .25°,  A,  =  1.Q5,  =  1.2275,  A^  =  0.9596. 


coefficient  and  a  as  a  function  of  time.  It  appears  that 
the  free  motion  a,  C^,  and  Cj^  are  decreasing  during  the  last 
two  cycles. 

Figure  26  shows  an  unstable  case  with  =  1.50° 
and  all  other  parameters  as  before.  It  is  observed  that  C 

m 

now  clearly  contains  a  higher  frequency  component.  The  lift 

time  history  of  Figure  27  does  not  show  a  similar  higher 

frequency  component.  These  results  are  typical  of  all  the 

cases  studied  by  the  time  history  method. 

In  Figure  2  8,  a  master  plot  is  shown  which  gives 

the  required  damping  parameter  needed  to  prevent  flutter  as 

a  function  of  the  imposed  disturbance  am.plitude.  Numerical 

integrations  are  indicated  by  the  solid  curve  at  intervals 

of  ct,  equal  to  .25°  up  to  1.5  0°.  A  Fourier  analysis  of  the 

last  forced  cycle  of  C  is  used  as  input  into  the  method  of 

m 

imposed  disturbances  to  generate  the  dashed  curve.  Both 

methods  are  shown  to  predict  the  same  general  trends.  The 

initial  effect  of  amplitude  on  required  damping  is  one  of 

a  moderate  decrease  up  to  about  =  .25°.  After  this,  system 

damping  must  be  continually  increased  in  order  to  prevent 

flutter  at  increasing  imposed  disturbance  amplitudes. 

Eleven  representative  points  ,  .  .  .  have 

been  selected  for  detailed  study  herein.  The  time  histories 

of  Figures  24  and  25  correspond  to  Point  4  and  Figures  26  and 

27  are  for  Point  6.  The  Fourier  coefficients  a.  a.  ...,a,,. 
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Figure  26a.  Moment  and  Angle  of  Attack  Time  History  for  NACA  64A006  Airfoil, 
M  =  .88,  a,  =  1.50°,  A,  =  1.05,  A_  =  1.2275,  A^  =  0.9596. 


Figure  28. 


ct^  (degrees) 

Required  Damping-Amplitude  Flutter  Boundary  for 
NACA  64A006  Airfoil  Pitching  in  Transonic 
Flow. 
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are  shown  in  Table  2  for  Point  4.  By  far,  the  a^ coefficient 
is  the  largest  for  both  and  C. .  This  indicates  that  the 

m  ij 

equation  of  motion  is  nearly  linear  at  a  disturbance  level 
of  =  .25°.  Table  3  shows  the  Fourier  coefficients  at 
=  1.50°  (Point  6) .  The  greatest  order  of  magnitude 
increase  between  Tables  2  and  3  is  in  the  a^  term  which  has 
increased  by  a  factor  of  about  250. 

The  results  of  a  sequence  of  Fourier  analyses 
is  given  for  Points  6  and  9  in  Tables  4  and  5,  respectively. 
The  parameter  being  varied  is  which  is  the  iime  after  which 
the  airflow  is  assumed  to  be  periodic.  A  Fourier  analysis 
of  the  first  cycle  of  imposed  motion  ^/2i\=Q)  for  Point  6 
is  shown  to  underestimate  a^^  by  about  seven  percent, 
compared  to  a  cycle  4  analysis  (y'^/2Tr=3).  By  Equation  (4.35) 
this  also  leads  to  a  value  of  required  damping  which  is 
fifteen  percent  low.  For  Point  9,  the  trend  is  reversed 
and  a  first  cycle  analysis  overestimates  a^  and  hence  by 
about  13%. 

The  phase  plane  plot  has  shown  to  be  a  useful 
way  to  determine  stability  for  this  problem.  As  can  be  seen 
from  the  time  history  plots  in  Figure  24,  it  can  be  difficult 
sometimes  to  judge  whether  or  not  the  amplitude  of  the  free 
motion  is  decaying.  Figure  29  shows  the  (ct,a')  phase  plane 
plot  for  Point  4.  The  circle  of  radius  .25  is  generated 
during  the  imposed  sinusoidal  motion.  The  free  motion  starts 


TABLE  2 


FOURIER  COEFFICIENTS  FOR  C. ,C 

L  in 

POINT  4.  Y,  =  6n 


a. (C^ ) xlO^ 

X  Li 

a. (C„) xlO^ 

1  m 

.70929 

-.04751 

-5.61950 

.44031 

-  .00006 

.00052 

.00317 

-.00438 

-  .00022 

-.00017 

.00070 

.0004^ 

-  .00011 

.00009 

.00064 

.00009 

.00032 

-.00003 

-  .00011 

.00017 

.00032 

.00019 

TABLE  3 


[>X*^ 


TABLE  4 


CONVERGENCE  OF  a,  FOR  C 

1  m 

POINT  6 


a,  (C„) xl0‘ 
1  m 


0 

1 

2 

3 


3.06340 

3.24630 

3.29070 

3.30140 


TABLE  5 

CONVERGENCE  OF  a,  for  C 

1  m 

POINT  9 


Yj^/2tt 

a, (C„) xlO^ 

1  m 

0 

2.24520 

1 

2.04880 

2 

2.00330 

3 

1.99330 

100 


ALFONE  -  .  250 

A1  -  1.  0500 

A2  -  1.2275 

A3  -  .  9595 

Figure  29.  Phase  Plane  Plot  for  Point  4  on  Dcimping- 
Amplitude  Plot. 
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at  the  end  of  the  fourth  imposed  cycle  ata=0,  a'  =  ,25 
("START"  on  the  figure).  It  can  be  seen  immediately  after 
the  free  motion  starts  that  the  amplitude  is  decreasing. 

Such  a  trend  is  more  difficult  to  observe  on  the  time 
history  plot.  Figure  30  shows  phase  plane  for  the  unstable 
Point  9 .  Phase  plane  plots  for  the  other  represe.itative 
points  on  Figure  28  are  included  as  Appendix  A. 

Representative  pressure  coefficient  distribu¬ 
tions  generated  by  LTRAN2  are  shown  in  Figures  31  through 
37.  The  distributions  are  for  the  last  cycle  of  imposed 
motion  for  Point  4  on  Figure  28.  At  a  =  45°,  Figure  30 
shows  the  existence  of  a  region  of  local  supersonic  flow 
between  10%  and  70%  of  the  chord  on  the  upper  surface. 

The  region  is  indicated  by  the  fact  that  Cp^  exceeds  the 
value  of  the  critical  pressure  coefficient,  Cp*,  in  this 
region.  A  shock  is  located  approximately  at  the  point  where 
the  nearly  discontinuous  drop  in  Cp  to  values  lower  than 
Cp*  occurs.  In  Figure  31,  the  shock  on  the  lower  surface  is 
slightly  to  the  rear  of  the  one  on  the  upper  surface.  At 
a  =90°  as  shown  in  Figure  32,  the  upper  surface  shock  has 
moved  aftward  and  the  lower  surface  shock  has  moved  forward. 
This  trend  continues  in  Figure  33  at  a  =  135°  where  it  is 
noted  that  the  lower  surface  shock  starts  to  weaken.  Little 
change  can  be  seen  in  Figure  34  at  a  =  180°;  however,  at 
a  =  225°  Figure  35  shows  the  lower  surface  shock  is  now 
moving  aftward  and  the  upper  one  is  moving  forward.  The 


steady  Pressure  Coefficient  Distribution,  NACA  64A006 
rfoil  Pitching  in  Transonic  Flow,  =  .88,  n  =  .25° 


Figure  35.  Unsteady  Pressure  Coefficient  Distribution,  NACA  64A006 
Airfoil  Pitching  in  Transonic  Flow,  =  .88,  a  =  .25° 
a  =  225°.  ”  ^ 


Figure  37.  Unsteady  Pressure  Coefficient  Distribution,  NACA  64A006 
Airfoil  Pitching  in  Transonic  Flow,  M  =  .88,  a  =  .25° 


Figure  38.  Unsteady  Pressure  Coefficient  Distribution,  NACA  64A006 
Airfoil  Pitching  in  Transonic  Flow,  M  =  .88,  a  =  .25“ 
n  =  360“.  “  ^ 


shock  positions  now  pass  each  other  at  about  270°  into  the 
cycle  as  shown  in  Figure  3G .  The  aftward  motion  of  the 
Cpj^  shock  and  forward  motion  of  the  Cp^  shock  continues  in 
Figures  37  and  38  which  concludes  one  cycle.  These  figures 
show  that  the  largest  change  in  the  pressure  distribution 
occurs  in  the  region  of  the  shock  which  has  a  mean  position 
of  about  75%  chord. 

4.4  TWO-DEGREE-OF-FREEDOM  FOLLOWER  LOAD  STUDY 

The  stability  of  a  bar  restrained  at  one  end  and 
compressed  at  the  other  end  by  a  "follower"  force  has  long 
been  a  problem  of  interest  in  analytical  mechanics.  A 
follower  force  is  a  compressive  load  which  after  deformation 
rotates  together  with  the  end  section  of  the  bar.  At  all 
times  it  remains  tangent  to  the  deformed  axis.  The  book  by 
Bolotin  gives  a  quite  detailed  analysis  of  this  problem 
from  the  standpoint  of  a  continuous  bar  model.  The  object¬ 
ive  of  the  follower  load  problem  is  to  determine  the  maximum 
load  value  which  the  bar  can  sustain  without  unbounded  ampli¬ 
tude  motion  resulting.  This  is  a  fundamental  example  of  a 
nonconservative  elastic  stability  problem.  It  is  termed 
nonconservative  since  the  position  dependent  loading  cannot 
be  derived  from  a  potential  function. 

In  what  follows,  a  discrete  two-degree-of-freedom 
follower  load  problem  is  analyzed.  A  nonlinearity  is  added 
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to  the  problem  by  requiring  that  the  load  not  only  "follow" 

but  that  a  component  of  its  magnitude  be  dependent  upon  the 

square  of  the  slope  at  the  tip.  Such  a  nonlinear  stability 

problem  is  not  in  such  standard  works  on  nonconservative 

systems  as  the  Bolotin  reference.  A  literature  search  has 

8  3 

shown  that  Roorda  and  Nemat-Nasser  are  the  first  to 
present  an  energy  based  method  for  this  nonlinear,  non¬ 
conservative  problem.  The  Reference  83  linearized  solution 
is  used  below  in  an  application  of  the  method  of  imposed 
disturbances,  and  excellent  agreement  with  the  Roorda 
and  Nemat-Nasser  nonlinear  load-amplitude  curve  is  obtained. 
4.4.1  Derivation  of  the  Equations  of  Motion 

The  two-degree-of-f reedom  follower  load 
problem  is  shown  in  Figure  39.  The  bar  segments  A-B  and 
B-C  are  assumed  rigid  and  are  of  length  1.  Elastic  and 
damping  properties  of  the  bar  are  modeled  as  two  torsional 
springs  located  at  positions  A  and  B.  The  restoring 
moments  at  A  and  B  are,  respectively, 

(4.36a)  =  B^^j_ 

(4.36b)  Mg  =  c((J)2-!t)^)  +  B2(02~<^i) 
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where  c  is  an  elastic  spring  constant,  and  b2  are 
damping  constants,  and  and  $2  the  angular 

rotations  at  A  and  B,  respectively.  A  mass  m  is 
concentrated  at  B  and  a  mass  2m  at  C.  The  load  y  is 
always  directed  along  the  line  B-C  and  is  related  to 
the  tip  slope  by 

(4.37)  y  =  A  (l+a-f>2^) 

93 

Lagrange's  equations  (see  Meirovitch, 
p.  50,  for  example) 


(4.38) 


d  3L 


3L 

34)^ 


Q.NC 

1 


i  =  1,2 


are  used  to  derive  the  equations  of  motion  where  L  =  T-V 

NC 

is  the  Lagrangian  and  are  the  nonconservative  general¬ 

ized  forces.  The  kinetic  energy,  T,  can  easily  be  written 


as 


(4.39)  T  =  j  [34)2^^+4)2^+2())j^4)2  cos  (0j^-(t)2)  1 

and  the  potential  energy  due  to  the  spring  constants,  V, 

is 


(4.40) 


V  =  ^(20^2 


-  20^02  +  <t>2  ) 
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The  generalized  nonconservative  forces  are  given  by 

NC  2 

(4.41)  Q,^  =  AJl(l+a<|)2  )  sin ( 

and 

(4.42)  =  b2(i;L"'^2^ 

In  writing  the  equation  of  motion  it  is  found 
convenient  to  introduce  the  following  nondimensional 
variables 

(4.43a)  Y  =  /c/m£^  t 

(4.43b)  =  b^/2,v/^  i  =  1,2 

(4.43c)  F  =  XJl/c 

Substitution  of  Equations  (4.39)  through 

(4.43)  into  Equations  (4.38)  results  in  the  nonlinear 
equations  of  motion 

2 

3(5j^  +  $2  cos((})j^-())2)  +  sin((J)^-<{)2) 

(4.44a)  ,  _ 

+  2()>^-((>2  +(Bj^+B2)  iJ2_“B2'i>2  ~  F(l+a({)2  )  sin ( 4)j^-<)!2) 

'<^2  'fj^cos  (<J)^-(J)2)  -  ())^+())2-B24>j_+B2<j'2 
(4.44b)  .  2 

sin(  (})j^-(j)2)  =  0 

where  the  "dot  notation"  for  derivatives  now  represents 
differentiation  with  respect  to  y •  If  all  the  nonlinear 
terms  in  Equations  (4.44)  are  set  to  zero,  the  linearized 
equations  can  be  written  in  matrix  fc’^m  as 
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3  I  N) 


1  IN 


(4.45) 


B1+B2 


2  -n  u 


B2  $2  ^  ^^2 


=  F 


1  -1 


0  0  u 


Following  the  notation  of  Reference  83 ,  a 
steady-state  solution  of  Equation  (4.45)  can  be  written 


(4.46) 


=  A  I  2^  cosay  -  2^ 


sin^Y 


where  A,  Z  ,  and  Z,  are  constants,  and  a  is  the  nondimen- 

ci  O 

sional  linearized  solution  frequency .  Equations  (4.46)  are 
substituted  into  Equations  (4.45)  and  the  coefficients  of 
cosay  and  sinfiy  set  to  zero  to  obtain  a  fourth  order 
set  of  homogeneous  equations  in  the  unknown  parameters  a 
and  F.  By  setting  the  determinant  of  the  system  equal  to 
zero,  a  characteristic  polynomial  in  a  results.  Refer¬ 
ence  83  solves  this  polynomial  and  obtains  the  following 
linearized  solution 


(4.47a) 

= 

0 

(Bj^+B2)  /(B^+6B2) 

(4.47b) 

Z  = 
a 

i-2a^^(a^^-i)/[(a^^-i) ^+a^^B2^] 

(4 .47c) 

^b  = 

3  2222 

-2aQ^B2/[(fio  -1)  ^+a^^B2^1 

(4  .4  7d) 

F  = 
0 

(4Bj_^+33B^B2+4B2^)  /2(Bj^+B2)  (B^+6B2)  +B^B2/2 
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derived  from  Hamilton's  principle.  Their  procedure  is 
equivalent  to  conserving  energy  in  the  limit  cycle  tto  qb-^in 
A)  and  requiring  that  the  equations  of  motion  be  satisfied 
(to  get  Z  ,  and  Z,  )  . 

Si  D 

The  method  of  imposed  disturbances  requires 
only  that  the  condition  =  0  is  satisfied  where  is 

a  function  of  some  appropriate  response  mode;  i.e.,  imposed 
disturbance.  Numerical  values  for  this  problem  can  be 
obtained  if  B^<<B2<<1,  products  of  B^  and  B2  are  small,  and 
if  3  E  B^/B^  is  given.  For  6  =  10,  Equation  (4.50)  is  used 
to  compute  F  =  F(A)  where  il,  Z  ,  Z,  are  the  linearized  values 

Si  D 

obtained  from  Equations  (4.47)  .  The  expression  for 
a  =  1/2  is 

F  =  24. 34/(11. 67+88. 61A^)  , 

and  for  a  =  0  it  is 

F  =  24 .34/(11. 67-28. 25A^) 

The  results  are  shown  in  Figure  4  0  and  compared 
to  those  of  Reference  33  .  For  ot  =  1/2  the  system  has  what 
was  termed  in  Chapter  3  as  "soft  nonlinearity."  A  "hard 
nonlinearity"  occurs  when  a  =  0.  The  results  of  both  methods 
agree  quite  well  especially  below  amplitudes  of  about  0.3. 

4.5  HYPERSONIC  PANEL  FLUTTER  STUDY 

An  elastic  panel  subjected  to  hypersonic  airflow  over 
one  side  can,  under  certain  combinations  of  flow  and 
structural  parameters,  undergo  large,  destabilizing 
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oscillations.  In  this  section,  an  elastic  panel  of  infinite 
span  is  studied  in  two-dimensional  nonlinear  hypersonic 
flow.  The  panel  is  mounted  on  a  viscoelastic  foundation. 
Boundary  conditions  are  simple-suuports  and  an  edge  compres¬ 
sion  less  than  the  static  buckling  load  is  allowed.  Only 
aerodynamic  nonlinearities  are  included.  The  Galerkin 
method  is  used  to  arrive  at  a  system  of  ordinary  nonlinear 
equations.  A  static  subcase  solution  of  the  method  of 
imposed  disturbances  is  compared  to  numerical  integration, 
and  it  is  shown  that  both  methods  predict  similar  load- 
amplitude  curves.  It  is  then  shown  how  the  method  of 
imposed  disturbances  can  be  used  to  perform  parametric 
studies.  The  parameter  selected  is  the  foundation  stiff¬ 
ness,  and  it  is  shown  that  great  computational  efficiency 
can  be  achieved  using  the  method  of  imposed  disturbances 
compared  to  multiple  time  history  solutions. 

4.5.1  Galerkin  Solution  of  Panel  Equation 

The  simply  supported  sandwich  panel  of  length 

I  is  shown  in  the  sketch  on  the  next  page.  A  hypersonic 

air  flow  is  loading  the  upper  surface  of  the  panel.  A 

compressive  force  is  applied  to  the  ends  and  the  lower 

surface  is  mounted  on  a  viscoelastic  foundation  The 

structural  properties  of  the  panel  are  the  bending  stiff- 

94 

ness,  D,  and  the  mass  density,  m.  The  book  by  Allen 
gives  expressions  for  D  in  terms  of  sandwich  geometric  and 
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material  properties.  It  will  be  sufficient  for  the 
purposes  herein  to  state  that  a  set  of  design  variables 
such  as  face  sheet  thickness,  core  thickness,  and 
elastic  material  constants  can  be  used  to  define  all  the 
structural  parameters.  The  results  presented  are  for  non- 
dimensional  groupings  of  constants,  and  if  needed,  these 
terms  can  be  expressed  as  functions  of  the  basic  sandwich 
panel  properties.  Such  a  procedure  is  necessary  to  solve 
a  specific  sandwich  panel  flutter  problem,  but  it  does 
not  add  further  insight  into  the  development  of  the 
method  of  imposed  disturbances. 

The  governing  nonlinear,  partial  differential 
equation  for  the  panel  is 


Dw  + 

xxxx 


Pw 


XX 


+  mw 


tt 


=  -|2[w  w 

X  V_  t  4 


(4.51) 


-  kj^w  -  k2W^ 


where  w  is  the  transverse  deflection  of  the  panel,  P  is 
the  compressive  edge  load,  q  is  the  dynamic  pressure  of  the 


airstream,  is  the  free-stream  Mach  number,  is  the  free- 
stream  velocity,  <  is  the  ratio  of  specific  heats  (1.4  for 
air)  ,  and  and  k2  are  the  elastic  and  damping  constants 
of  the  foundation.  The  subscripts  x  and  t  denote  partial 
differentiation  with  respect  to  the  streanwise  coordinate 
and  time,  respectively.  The  left  side  of  Equation  (4.51) 
is  derived  from  elementary  beam  theory,  and  the  left  side 

airloading  terms  come  from  a  modified  second  order  piston 
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theory  expansion  for  the  pressure.  McIntosh  states 

2 

that  the  nonlinear  terms  w  and  w  w  are  the  two  m.ost 

A  X  ^ 

important  for  hypersonic  flow. 

It  is  convenient  to  introduce  the  following 
nondimens ional  variables 

C  =  x/Z 

T  =  /^)/M  I*  t 

CO 

w  =  w/1 
"u  =  /£p  /mM 

^1  “ 

k2  =  / I'^/mD  k2 
P  =  l^P/D 
X=  V  ^/DM^ 

CO  CO  '  CO 

M  =  (<+l)M  /4 
00 
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so  Equation  (4.51)  becomes  in  nondimensional  form 


(4.52a) 


P  w^_+  k,w  +  )c_w  +  w 

1  2  T  TT 

=  -A^(w^  +M  ~  Xu(w^  +  2M  w^w^) 


on  0<5<1.  The  boundary  conditions  for  the  simply  supported 
condition  are 


(4.52b)  w  =  0 

(4.52c)  =  0 


at  ?  =  0,1 


The  Galerkin  solution  (see  Chapter  6  of 
Reference  93,  for  example)  has  the  form 


_  N 

(4.53)  w(C,t)  =  I  (T)sinjTT5 

j=l  ^ 

where  q^  are  as  yet  undetermined  functions.  The  terms 
si  ijTTC  are  seen  to  satisfy  all  the  boundary  conditions  for 
the  problem.  Equation  (4.53)  is  substituted  into 
Equation  (4.52a>,the  result  is  premultiplied  by  siniir^,  and 
an  integration  on  C  from  0  to  1  is  performed.  The  resulting 
N  equations  are 
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;4  4  -  ^2^2  -  ,  ,(1) 


]  IT  -P  j  TT  +k 


(4.54) 


+  I  Ilj  qj  =  I  (-A^jTT)!^])  +  I  (-A.^)l|l)q 

j=l  J  j=l  13  j  ID  D 

+  I  I  (-X^M  jk^^)  q  q 

j=l  k=l  ^  ^ 


N  N 


+  I  i  (-2Xu  MjTT)  q  q 

j=l  k=l  ijk  => 


The  Galerkin  integrals  and  are  defined 

in  Appendix  B.  It  is  shown  that  l|V  is  equal  to  .  j  where 

6.  .  is  the  Kronecker  delta.  This  expression  for  is 

ID  ID 

incorporated  into  the  first  order  equations  which  follow. 


(4.55a)  x^  =  q^ 


(4.55b) 


'i+N  *^i 


so  Equations  (4.54)  have  the  final  form 
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=-(k2+XTr)x.^j^ 


(tt  1  -Fit  i  +k^) 


i+N 


(4 .56) 


-2  A  TT 


N 


I  j  lf2)x. 
jil  ^ 

7-  7  N  N  , 

-2A  ^  I  j  k  X.  X 

j=l  k=l  ^ 

—  r  ^  (4) 

■  ’'ill  Ji  ^  ’‘3 


Equations  (4.56)  total  2N  in  number  and  can  be  solved  by 
numerical  integration. 

4.5.2  Calculations  of  Nonconservative  Work 


The  total  nonconservative  work  done  on  the 
panel  by  the  hypersonic  flow  and  the  viscoelastic  foundation 
during  one  period  of  steady-state  motion  is 


(4.57) 

= 

-AyJ^ 

where 

T 

1 

(4.58a) 

j 

r 

0  - 

w-w  dCdx 

0  ^ 

(4.58b) 

1 

[  w,^  w  d^dT 

> 

0  J 

1  0  T 

(4.58c) 

rT 

[  ^w  ^  dCdT 

0  ■ 

'0  ^ 

(4 .58d) 

•T 

{  ^  -  2 

w_w  dCdx 

Using  the  Galerkin  solution  Equation  (4.5  3)  ,  these  integrals 
can  be  rewritten  as 


(4.59a) 

(4.59b) 

(4 .59c) 

(4 .59d) 

where 

(4.60a) 

(4.60b) 

(4.60c) 

(4 .60d) 


.(1) 


.(2) 


.(3) 


.(4: 


N  N 

=  -  I  !  j 

itl  5=1  13  13 


N  N  N 


i=l  j=i  ic=l  ^3K 

N 


=  I  ^  ''i 
^  i=l  ^ 


(3) 


N  N  N  ... 

=  ^  I  I  I  j  ij jk 

i=l  j  =  l  k=l 


ID 

.(2) 

^ijk 


k!2)  = 


Kf3)  = 
1 


= 

illc 


T  . 

qi  q.  dT 
0  ^  3 

T 

q  q.  q  dx 
0  1  3 

T 


*  2  , 
q.  dx 


^i 


dx 


4.53  Application  of  the  Method  of  Imposed 
Disturbances 


Let  the  generalized  coordinates,  q^ ,  have  the 


form 

(4.61)  q^(x)  =  A4'^(x) 

where  A  is  an  undetermined  amplitude  pareimeter  and  t  is  an 
appropriate  imposed  disturbance.  Equation  (4.57)  can  now  be 
written 
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ii 


(4.61) 


NC  2 
-W^^/A^  = 

where 

(4.62a)  =  X^M  +  2X11  M 

(4.62b)  =  X^  +  Xy 

and  the  J  quantities  are  evaluated  using  the  K  integrals  of 

Equations  (4.60)  which  have  the  q's  replaced  by  fs. 

NC 

At  steady  state  motion,  W  =0,  and  X  =  A 
and  so  Equation  (4.61)  is  solved  to  get  A  in  terms  of  A 
and  the  system  variables  to  get 

(4.63)  A(A)  =  -T2(A)/T^(A) 

To  determine  the  function  explicitly,  it  is 
now  necessary  to  select  an  imposed  disturbance,  1'.  The  com¬ 
ponents  of  the  imposed  disburbance  are  selected  to  be  the 
steady-state  solution 

(4.64)  =  a^siny  i=l,N 

where  y  =  .  Substitution  of  these  quantities  into 

Equations  (4.60)  for  the  q  terms  gives 

(4.65a)  =  0 

13 

(4.65b)  =  0 

(4.65c)  i=l,N 

1  1 

(4.65d)  =  0 
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and  thus  by  Equations  (4.52) 


(4.66a)  =  0 

and 

-  -  N  , 

(4.66b)  T-  =  (7^tj/2)  (Xu+k_)  I  a. 

^  ^  i=l  1 

Substitution  of  and  into  Equation  (4.61)  gives  the 
result 

N 

(4.67)  =  -(uwA^/2) (Xu+k  )  I  a.^  <  0 

^  i=l  ^ 

This  implies  that  there  is  no  steady-state,  limit  cycle 
type  instability  based  on  the  selected  imposed  disturbance, 
Equation  (4.64) .  It  will  be  shown  in  the  next  subsection 
that  the  static  subcase  of  the  method  of  imposed  distur¬ 
bances  does  lead  to  an  amplitude  dependent  instability 
which  occurs  at  zero  frequency. 

4.5.4  The  Static  Subcase  of  Nonlinear  Panel 
Instability 

The  case  of  static  instability  is  now  checked. 
As  the  starting  point,  assume 
_  _  ^ 

(4.68)  w  =  w(C)  =  ^  q-  sinjTT^ 

3  =  1  ^ 

where  the  q's  are  constants.  Application  of  the  Galerkin 
method  as  shown  in  Section  4.5.1  leads  to  the  following 
system  of  nonlinear  algebraic  equations 
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(4.69) 


j=l  ^  -  3 


.2 
D  q 


N 


ifl) 

:  i: 

(2) 


•  -  I 
iii 

+  ic,  y  q  .  +  X'^TT  y  j  q. 

lj=l  ^  j=l  1  13 

P-  7  N  N 

+  X^M  TT^  y  I  j  k  q  q  =  0 

j=l  k=l  ^ 

For  N=l,  these  nonlinear  force  balance 
equations  reduce  to 


(4.70a) 


ZF  =-f„  +  f  =  0 

S  3. 


where  the  restoring  force,  f^,  is 


(4.70b) 


,  1,4  2-  ,  r-  , 

fg  =  -  TT  P  +  k^) 


and  the  applied  load  due  to  the  flow  is 
(4 .70c)  f,  =  -I  X^  Mirq,^ 

Solving  equations  (4.70)  for  q^^  gives 
(4.71)  q^  =  -3(Tr'*-Tr^P  +  k^)/4X^M  tt 

This  last  equation  is  an  approximation  to  the  required  load- 

_  2 

amplitude  plot.  It  is  seen  that  for  P<tt  (the  critical 

Euler  buckling  load) ,  q^  will  be  negative.  Thus  for  any  other 
*  * 

q^,  e.g.,  q^  such  that  q^<q^<0,  there  will  be  a  force 

imbalance  in  Equation  (.4.70a)  and  the  panel  will  collapse 

6  9 

inward.  Such  a  failure  mode  is  noted  by  McIntosh. 
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4.5.5  Numerical  Results 


The  results  of  a  computer  analysis  of  the 
hypersonic  panel  are  now  presented.  The  code  performs  a 
Kutta-Merson  numerical  integration  of  the  equations  of 
motion.  Equations  (4.56),  and  has  the  following  graphical 
output  displays  available: 

(1)  Phase  plane  plots  for  any  point  on  the 
panel  (w^  vs.  w  at  a  specified  point) 

(2)  Time  history  of  displacement  and  velocity 
at  a  given  point  (w,  w^  vs.  t  at  a  point) 

(3)  Mode  shape  plots  at  a  given  time  (w  vs.  C 
at  a  given  t) 

(4)  Time  history  of  the  total  work  done  by  the 

external  forces  acting  on  the  panel.  This  quantity,  which 

is  denoted  by  WTD  on  the  plots,  can  be  integrated  over  one 

cycle  of  motion  to  get  the  total  nonconservative  work, 

NC 

W  ,  given  by  Equation  (4.5  7)  . 

The  code  allows  for  up  to  10  modes  (N=10)  in 

the  Galerkin  solution. 

A  separate  code  solves  the  nonlinear  force 

balance  equations  (Equations  4.69))  by  a  standard  Newton- 
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Raphson  method  (see  Ralston,  Chapter  8,  for  example)  to 
get  a  set  of  q's.  These  particular  q's  are  referred  to 
as  and  can  be  written  as 

(4.72) 


L3] 


where  A  is  the  maximum  and  4'^  are  normalized  constants. 
The  ^ ' s  are  input  into  the  Kutta-Merson  integration  program 
as  initial  displacement  conditions.  By  selecting  various 
A  values  at  the  given  X  values,  a  stability  boundary  can  be 
determined  which  gives  the  X  =  A (A)  curve.  It  is  shown 
that  the  value  of  A  predicted  by  the  Newton-Raphson  method 
agrees  quite  well  with  the  one  obtained  parametrically  in 
the  Kutta-Merson  code. 

Figures  41  through  44  show  typical  stable 
panel  motions  (Case  1) .  Data  for  this  case  is  given  on 
the  phase  plane  plot  of  Figure  41.  In  Figure  42,  the  time 
histories  of  w  and  w^  are  shown.  The  entire  panel  mode 
shape  at  T  =  0  is  given  in  Figure  43  which  is  followed  by 
the  work  input  time  history.  These  plots,  taken  together, 
show  that  the  panel  is  stable  at  an  amplitude  level  of  -.04. 

Typical  stable  motion  for  N=2  is  shown  in 
Figures  45  and  46.  At  an  amplitude  level  of  -.04  the  motion 
is  stable  as  evidenced  by  the  phase  plane  plot  (Figure  43) 
and  work  input  time  history  (Figure  46) .  By  just  changing 
the  amplitude  to  -.0588,  the  unstable  Case  3  phase  plane 
plot  is  obtained  in  Figure  47.  The  work  input  is  clearly 
shown  to  be  increasing  in  the  time  history  of  Figure  48 . 


^1  t:| 


CASE  L: 
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50 


HYPERSONIC  PANEL  FLUTTER 
X/L  -  .  750 

H  -  .010  NSTEP  -  150 

CPTIME  -  .  517 


Figure  41.  Phase  Plane  Plot  Showing  Stable  Panel 
Motion,  Case  1. 


HYPERSONIC  PANEL  FLUTTER 
X/L  =  .  750 

H  -  .010  N3TEP  -  150 


Figure  46.  Time  History  of  Work  Input  to  Panel 


Table  6  shows  a  systematic  variation  of  all 
the  parameters.  The  static  subcase  of  the  method  of 
imposed  disturbances  is  used  to  generate  4'^,  and  the 

amplitude  The  normalized  terms  4'^,  4*2  are  treated  as 

imposed  displacements  in  the  time  history  solution.  These 
conditions  are  then  scaled  until  the  motion  is  neutrally 
stable .  The  scale  factor  computed  by  Kutta-Merson  inte¬ 
gration  is  and  is  shown  in  the  last  column  of  Table  6. 

The  effect  of  edge  compression  is  seen  to  lower  the  ampli¬ 
tude  at  which  the  system  first  becomes  unstable.  On  the 
other  hand,  increasing  the  stiffness  of  the  foundation  has 
the  beneficial  effect  of  increasing  the  amplitude.  Foun¬ 
dation  damping  is  shown  to  have  little  effect  on  amplitude. 
Increasing  the  Mach  number  parameter,  M,  lowers  the  critical 
amplitude  but  the  mass  ratio,  u,  variation  has  no  effect. 

In  all  cases  shown  in  Table  6 ,  the  method  of  improved 
disturbances  and  Kutta-Merson  integration  give  very  compar¬ 
able  results. 

The  relative  computing  efficiency  of  these 
methods  is  shown  in  Tcibles  7  and  8.  Table  6  shows  the 
normalized  imposed  disturbances,  the  amplitude  parameter, 
and  the  computing  time  needed  to  solve  the  N=4  static 
equation  by  the  Newton-Raphson  procedure.  Table  7  shows 
results  for  numerical  integration  using  the  predicted 
disturbances  and  two  selected  amplitudes,  and  A2 .  The 
first  is  slightly  lower  than  A^^  and  the  second  slightly 
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TABLE  7 

METHOD  OF  IMPOSED  DISTUPEANCES 
FOUNDATION  STIFFNESS  STUDY 
N=4,  M=5,  11’=.  2,  P=k-=0,  k,  =40 


'^2 

^3 

^ID 

-.04661 

.02855 

-.00275 

-.29213 

-.06722 

.03970 

-.00401 

-.20583 

-.12030 

.03443 

-.00748 

-.12235 

-.19054 

.04410 

-.01274 

-.08598 

-.28087 

.06087 

-.02072 

-.06693 

-.39517 

.08705 

-.03233 

-.05427 

-.53671 

.12576 

-.04821 

-.04387 

-.70736 

.18192 

-.06919 

-.03482 

TABLE 

8 

KUTTA-MERSON  NUMERICAL  INTEGRATION 

FOUNDATION  STIFFNESS  STUDY 

N=4,  M=5 

,  y=.2, 

P=k2=0,  k 

1=40 

^  1  A, 

1  A, 

1  Comment 

iCP  Time 
1  / _ \ 

-.293 
-.207 
-.123 
-.087 
-.068 
-.055 
-.044 
-.0  36 


CP  Time 
( sec) 


.0  39 


62.889 

62.594 

62.548 

58.121 

61.056 

53.908 

52.239 

50.650 


higher.  The  conmient  column  of  Table  8  shows  that  the 
lower  A  produces  stable  motion  (S)  and  the  higher  gives 
unstable  motion  (U) .  Thus  the  instability  is  bracketed 
in  amplitude.  Comparing  the  computing  times,  it  can  be 
seen  that  for  this  problem  the  method  of  imposed 
disturbances  is  almost  2000  times  faster  than  numerical 
integration  and  has  equivalent  accuracy  in  predicting 
load-amplitude  behavior. 

Figure  49  shows  the  load-amplitude  plots 
N=l,2,3  and  4.  In  all  cases  the  method  of  imposed 
disturbances  and  numerical  integration  give  good  agreement 
for  amplitudes  above  .05.  Below  this  amplitude,  the 
method  of  imposed  disturbances  should  be  supplemented  by 
a  linear  flutter  analysis  which  predicts  the  N=2 ,  low 
amplitude  flutter  very  accurately. 


N=4 ,  Method  of  Imposed 
Disturbances 

N=2,  Method  of  Imposed 
Disturbances 

N=l,  Method  of  Imposed 
Disturbances 
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SECTION  V 


CONCLUSIONS  AND  RECOMMENDATIONS 

This  research  effort  has  resulted  in  an  efficient 
and  easy  to  implement  procedure  for  the  stability  assessment 
of  nonlinear,  nonconservative  systems.  The  importance  of 
the  analytical  procedure  is  discussed  in  the  first  section 
below.  The  chapter  then  concludes  with  a  section  describ¬ 
ing  areas  where  additional  work  is  needed.  Further  method 
development  areas  and  applications  are  indicated  in  that 
section. 

5.1  CONCLUSIONS 

The  stability  assessment  procedure  developed  herein 
is  based  on  the  principle  of  conservation  of  energy  in  the 
limit  cycle.  By  imposing  a  motion  on  a  system  and  treating 
the  amplitude  of  motion  as  a  variable,  the  total  non¬ 
conservative  work  done  by  the  nonlinear  position  dependent 
external  forces  and  the  internal  damping  forces  is  easily 
computed.  From  the  condition  that  the  total  nonconserva¬ 
tive  work  is  zero  during  steady-state  nonlinear  motion,  a 
relationship  among  the  dynamic  load  parameter,  the  amplitude 
of  motion,  and  the  system  variables  can  be  established. 


146 


Such  information  is  required  in  order  to  perform  optimiza¬ 
tion  studies  on  nonlinear  systems. 

The  method  of  imposed  disturbances  has  a  clear 
advantage  over  current  analytical  methods  for  the 
following  reasons: 

-  a  definite  answer  to  the  stability  question  is 
obtained, 

-  numerical  instability  problems  associated  with 
time  history  solutions  are  avoided, 

-  an  explicit  analytical  form  for  the  external 
forces  is  not  required, 

-  nonlinear  transonic  flutter  studies  can  be 
performed  without  coupling  the  aerodynamic 
calculations  to  the  structural  equations, 

-  the  method  is  applicable  to  nonlinear  external 
loadings  containing  both  circulatory  and 
dissipative  forces. 

These  positive  aspects  of  the  method  plus  its  computational 
efficiency  make  it  a  prime  candidate  for  inclusion  in  an 
automated  design  procedure  where  a  structure  or  structural 
component  is  subjected  to  nonlinear,  nonconservative  forces. 
The  procedure  developed  herein  is  a  key  analytic  tool 
which  makes  the  optimization  of  nonlinear  systems  feasible. 

A  static  subcase  of  the  method  of  imposed  disturb¬ 
ances  is  applicable  when  the  external  forces  contain  terms 


147 


which  are  independent  of  velocity.  By  solving  nonlinear 
static  force  balance  equations,  a  relationship  between 
the  load  parameter  and  the  amplitude  of  the  deformation 
can  be  established.  The  solution  of  the  nonlinear  static 
equations  is  much  more  economically  obtained  than  the  time 
history  results  when  the  prime  mode  of  instability  is 
static  in  nature. 

5 . 2  RECOMMENDATIONS 

Further  work  is  needed  to  demonstrate  the 
applicability  of  the  method  to  lifting  surface  flutter 
problems.  The  method  is  shown  herein  to  give  good  quali¬ 
tative  results  for  transonic  pitch  flutter,  but  applica¬ 
tions  are  now  needed  for  multiple  degree-of-f reedom  airfoils 
in  transonic  flow.  The  effect  of  structural  nonlinearities 
now  needs  also  to  be  included  in  the  analysis.  Such 
applications  present  no  conceptual  difficulty  for  the  method 
of  imposed  disturbances,  and  this  work  can  be  accomplished 
in  a  straightforward  manner. 

The  hypersonic  panel  flutter  application  needs  to  be 
studied  further.  Alternate  imposed  disturbances  should  be 
used  to  calculate  the  nonconservative  work  to  determine 
under  what  conditions  oscillations  with  finite  periods  can 
be  sustained.  The  effect  of  large  panel  deflections  should 
also  be  included.  In  all  cases  presented  the  edge  compres¬ 
sion  is  kept  below  its  Euler  buckling  value,  and  further 
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studies  are  now  needed  at  higher  edge  loadings  to 
determine  the  relationship  between  the  static  and  dynamic 
instability  phenomena. 

The  method  needs  further  development  to  establish 
a  relationship  between  the  initial  conditions  and  the 
critical  amplitude  of  motion.  In  this  way,  critical  initial 
velocities  and  displacements  can  be  included  in  the  design 
process.  Such  considerations  are  seen  to  be  important  in 
order  to  consider  gust  loadings  and  other  transient  aero¬ 
dynamic  phenomena. 

Theoretical  work  is  needed  to  develop  criteria  for 
the  selection  of  imposed  disturbances.  In  the  problems 
reported  herein,  satisfactory  results  are  obtained  using 
solutions  to  the  linearized  system.  The  method  can  also  be 
applied  using  steady-state  solutions  to  the  nonlinear 
problem.  It  appears  feasible  to  develop  an  imposed  distur¬ 
bance  method  using  multiple  disturbances  to  define  the 
load-amplitude  plot  in  a  piecewise  manner.  Time  history 
solutions  at  selected  amplitudes  should  be  performed  to 
aid  in  the  selection  of  imposed  disturbances. 

From  a  design  standpoint,  the  key  results  of  this 
study  are  that  the  addition  of  system  damping  is  more 
important  than  stiffness  or  mass  changes  in  alleviating 
limit  cycle  instability.  For  systems  which  exhibit 
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nonlinear  static  instabilities,  stiffness  variables 
are  of  prime  importance.  These  trends  need  to  be 
verified  further  by  additional  applications  to  airfoil 
flutter  and  divergence  problems. 
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APPENDIX  A 


TRANSONIC  PITCH  FLUTTER  PHASE  PLANE  PLOTS 

Additional  phase  plane  plots  are  shown  in  the 
Figures  50  through  58.  The  point  number  given  on  each 
figure  refers  to  the  damping-amplitude  plot  of  Figure  28. 
In  most  cases  such  as  Figures  50,  52,  54,  55,  57,  and  58, 
the  phase  plane  trajectories  are  either  converging  toward 
the  origin  or  diverging  away  from  it.  In  these  situations 
the  points  are  not  too  closely  situated  to  the  stability 
boundary  of  Figure  28.  In  Figures  51,  53,  and  57,  the 
motions  have  a  dual  converging  and  diverging  character. 
These  points  are  generally  close  to  the  stability 
boundary . 
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Figure  52.  Phase  Plane  Plot  for  Point  3  on 
Damping- Amplitude  Plot. 


164 


165 


.  9596 


Figure  54.  Phase  Plane  Plot  for  Point  7  on 
Damping-Amplitude  Plot. 
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Figure  55.  Phase  Plane  Plot  for  Point  8  on 
Damping-Amplitude  Plot. 
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Figure  56.  Phase  Plane  Plot  for  Point  9  on 
Damping-Amplitude  Plot. 
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Figure  57.  Phase  Plane  Plot  for  Point  10  on 
Damping-Amplitude  Plot. 
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APPENDIX  B 


GALERKIN  INTEGRALS  FOR  HYPERSONIC  PANEL 


The  constants  in  Equations  (4.54)  which  result  from 
the  integration  of  trigonometric  functions  are  defined 
as  follows. 


(1) 

(B.l)  I..  =  siniTT^sinjTrCdC 


(B.2)  ifj^  = 


sinirr^cos  jirCdC 


( 3)  ( 

(B.3)  i'-v  =  siniiT^cos^TT^coskTr^d^ 


= 


siniiT^cos  jTT^sinkTT^d^ 


By  use  of  elementary  trigonometric  identities,  these  can 


(  ! 


he  written  as 

(1)  i  I  ^  ^ 

(B.5)  l|V  =< 

L  0  if  i  j 

(B.6)  l[j^  =  [S(i+j)  +S(i-j)]/2TT 

(B.7)  I^^^  =  [S(i+j+k)  +  S  (i+j-k) 


+  S(i-j+k)  +  S(i-j-k)  ]/4tt 


(B.8)  ifjj^  =  [S(i+j-k)  -  S(i+j+k) 

+  S(i-j-k)  -  S(i-j+k)  ] /4Tr 
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